!c$$$c%-----------------------------------------------------------%
!c                                                                |
!C     Perturbed Normal Mode Analysis                             |
!c        using Elastic Network Potential                         |
!c        for Calclating Dynamics of                              |
!c        of PROTEIN/LIGAND Complex                               |
!c        with Alternative Binding Sites                          |
!c                                                                |
!c       By  Dengming Ming                                        |
!c       Los Alamos National Lab                                  |
!c       Fudan University                                         |
!c       July 2004                                                |
!c       Jan  2005  Revised by Ming                               |
!c       May  2005  add Backbone Enhance by Ming                  | 
!c       Feb  2006  developed Perturbation Calculation by Ming    |
!c       Nov  2013  add multiple chain interactions by Ming       |
!c                                                                | 
!c%---------------------------------------------------------------%
module dpa
  implicit none
  real*8,dimension(:),allocatable        :: ev0,vec0
  real*8,parameter    :: RPRECI=2.22045D-16, RBIGST=4.49423D+307
  integer,parameter    :: IZERO=0, IONE=1, ITWO=2, ITHREE=3,IFOUR=4
  real*8,parameter     :: ZERO=0.D0, ONE=1.D0, TWO=2.D0,HALF=0.5D0
  
  interface cdiag
     module procedure charmmdiag9
!     module procedure charmdiagq
  end interface cdiag
contains
  subroutine charmmdiag9(N,nmode,hessian,ev,vec)  ! 1-D allocatable hess, vec--optional
    implicit none
    integer,intent(in)    ::  n,nmode      !N IS THE DIMENSION OF DOF
    real*8,allocatable,dimension(:),intent(inout)  :: hessian, ev
    real*8,allocatable,dimension(:),optional,intent(inout)  ::  vec
    real*8,allocatable,dimension(:)  ::  vec2  !working dimension
    REAL*8  A(N+1),B(N+1),P(N+1),W(N+1),TA(N+1),TB(N+1),Y(N+1)
    integer               ::  NDD,NADD,ISALLOCATED
    !
    NDD=1                     !NDD is the first mode to be found. 
    NADD=NDD-1                
    !
    if(allocated(ev)) deallocate(ev)
    allocate(ev(n),stat=isallocated)          !create atom link for a residue
    if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR EV(N) IN DIAG');
    !
    if(present(vec)) then
       if(allocated(vec)) deallocate(vec)
       allocate(vec(n*nmode),stat=isallocated)     
       if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR VEC(N*NMODE) IN DIAG');
       call CHARMM_DIAGQ(N,NMODE,A,B,P,W,TA,TB,Y,NADD,hessian,ev,vec)
    else
       allocate(vec2(n*nmode),stat=isallocated)         
       if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR VEC(N*NMODE) IN DIAG');
       call CHARMM_DIAGQ(N,NMODE,A,B,P,W,TA,TB,Y,NADD,hessian,ev,vec2)
       deallocate(vec2)
    end if
    !
    return
  end subroutine charmmdiag9

  SUBROUTINE CHARMM_DIAGQ(NX,NFRQX,A,B,P,W,TA,TB,Y,NADD,DD,EV,VEC)
    !
    !     THIS ROUTINE IS A CONGLOMERATION OF GIVEN, HOUSEC, AND EIGEN
    !     WHERE THE BEST FEATURES OF EACH WERE KEPT AND SEVERAL OTHER
    !     MODIFICATIONS HAVE BEEN MADE TO INCREASE EFFICIENCY AND ACCURACY.
    !
    !   By Bernard R. Brooks   1981
    !
    !   NX      - ORDER OF MATRIX
    !   NFRQX   - NUMBER OF ROOTS DESIRED
    !   DD      - SECOND DERIVATIVE MATRIX IN UPPER TRIANGULAR FORM
    !   VEC     - EIGENVECTORS RETURNED (NX,NFRQX)
    !   EV      - EIGENVALUES RETURNED (NX)
    !   A,B,P,W,TA,TB,Y - ALL SCRATCH VECTORS (NX+1)
    !   NADD    - NUMBER OF LOWEST ROOTS TO SKIP
    !
    !
    implicit none
      !
      INTEGER NX,NFRQX,NAD
      REAL*8,INTENT(INOUT)   ::  DD(*)
      REAL*8,INTENT(INOUT)   ::  EV(*),VEC(*)
      REAL*8 A(NX),B(NX),P(NX),W(NX),TA(NX),TB(NX),Y(NX)
      ! 
      REAL*8 ETA,THETA,DEL1,DELTA,SMALL,DELBIG,THETA1,TOLER,ETAR
      REAL*8 RPOWER,RPOW1,RAND1,FACTOR,ANORM,U,ANORMR
      REAL*8 SUM1,BX,S,SGN,TEMP,XKAP,EXPR,ALIMIT,ROOTL,ROOTX,TRIAL,F0
      REAL*8 AROOT,ELIM1,ELIM2,T,EPR,XNORM,XNORM1,EVDIFF
      INTEGER N,NADD,NEV,NEVADD,NTOT,I,IPT,J,IJ,NN,MI,MI1,JI,JI2,II
      INTEGER ML,ML1,L,M,K,MJ,MJ1,NOMTCH,NOM,IA,ITER
      INTEGER J1,MK,MK1,KK
      !
      REAL*8 ANUMX
      !
      ANUMX=ZERO
      DO I=1,NX
         A(I)=ANUMX
         B(I)=ANUMX
         P(I)=ANUMX 
         W(I)=ANUMX
         EV(I)=ANUMX
         TA(I)=ANUMX
         TB(I)=ANUMX
         Y(I)=ANUMX
      END DO
      !
      ETA=RPRECI
      THETA=RBIGST
      !
      N=NX
      NEV=NFRQX
      NEVADD=NEV+NADD
      !
      DEL1=ETA/100.0
      DELTA=ETA**2*100.0
      SMALL=ETA**2/100.0
      DELBIG=THETA*DELTA/1000.0
      THETA1=1000.0/THETA
      TOLER=100.0*ETA
      ETAR=1.0/ETA
      RPOWER=8388608.0
      RPOW1=RPOWER*0.50
      RAND1=RPOWER-3.0
      !
      ! Find largest element.
      FACTOR=ZERO
      NTOT=(N*(N+1))/2
      DO I=1,NTOT
         FACTOR=MAX(FACTOR,ABS(DD(I)))
      ENDDO
      !
      ! Check for zero matrix.
      IF(FACTOR.LE.THETA1) THEN
         WRITE(*,811)
811      FORMAT(' WARNING FROM <DIAGQ>. Zero matrix passed.',&
              ' Identity matrix returned.')
         DO I=1,NEV
            EV(I)=ZERO
            IPT=(I-1)*N
            DO J=1,N
               IPT=IPT+1
               VEC(IPT)=ZERO
               IF(I+NADD.EQ.J) VEC(IPT)=ONE
            ENDDO
         ENDDO
         RETURN
      ENDIF
      !
      ! Compute norm of matrix
      FACTOR=ONE/FACTOR
      IJ=0
      ANORM=ZERO
      DO I=1,N
         DO J=I,N
            IJ=IJ+1
            U=(DD(IJ)*FACTOR)**2
            IF(I.EQ.J) U=U*HALF
            ANORM=ANORM+U
         ENDDO
      ENDDO
      !
      ! Scale the matrix
      ANORM=SQRT(ANORM+ANORM)/FACTOR
      ANORMR=ONE/ANORM
      DO I=1,NTOT
         DD(I)=DD(I)*ANORMR
      ENDDO
      !
      NN=N-1
      MI=0
      MI1=N-1
      !
      ! Perform trigiagonalization
      DO I=1,NN
         SUM1=ZERO
         B(I)=ZERO
         JI=I+1
         IPT=MI+I
         A(I)=DD(IPT)
         IPT=IPT+1
         BX=DD(IPT)
         JI2=JI+1
         DO J=JI2,N
            IPT=IPT+1
            SUM1=SUM1+DD(IPT)*DD(IPT)
         ENDDO
         IF(SUM1.LT.SMALL) THEN
            B(I)=BX
            DD(MI+JI)=ZERO
         ELSE
            S=SQRT(SUM1+BX**2)
            SGN=SIGN(ONE,BX)
            TEMP=ABS(BX)
            W(JI)=SQRT(HALF*(ONE+(TEMP/S)))
            IPT=MI+JI
            DD(IPT)=W(JI)
            II=I+2
            IF(II.LE.N) THEN
               TEMP=SGN/(TWO*W(JI)*S)
               DO J=II,N
                  IPT=IPT+1
                  W(J)=TEMP * DD(IPT)
                  DD(IPT)=W(J)
               ENDDO
            ENDIF
            B(I)=-SGN*S
            !
            DO J=JI,N
               P(J)=ZERO
            ENDDO
            ML=MI + MI1
            ML1=MI1-1
            DO L=JI,N
               IPT=ML+L
               DO M=L,N
                  BX=DD(IPT)
                  P(L)=P(L)+BX*W(M)
                  IF(L.NE.M) P(M)=P(M)+BX*W(L)
                  IPT=IPT+1
               ENDDO
               ML=ML +ML1
               ML1=ML1-1
            ENDDO
            !
            !
            XKAP=ZERO
            DO K=JI,N
               XKAP=XKAP+W(K)*P(K)
            ENDDO
            DO L=JI,N
               P(L)=P(L)-XKAP*W(L)
            ENDDO
            MJ=MI+MI1
            MJ1=MI1-1
            DO J=JI,N
               DO K=J,N
                  EXPR=(P(J)*W(K))+(P(K)*W(J))
                  DD(MJ+K)=DD(MJ+K)-EXPR-EXPR
               ENDDO
               MJ=MJ+MJ1
               MJ1=MJ1-1
            ENDDO
         ENDIF
         MI=MI+MI1
         MI1=MI1-1
      END DO
      !
      ! Begin sturm bisection method.
      !
      A(N)=DD(MI+N)
      B(N)=ZERO
      !
      ALIMIT=ONE
      DO I=1,N
         W(I)=B(I)
         B(I)=B(I)*B(I)
      ENDDO
      DO I=1,NEVADD
         EV(I)=ALIMIT
      ENDDO
      ROOTL=-ALIMIT
      !
      DO I=1,NEVADD
         ROOTX=ALIMIT
         DO J=I,NEVADD
            ROOTX=MIN(ROOTX,EV(J))
         ENDDO
         EV(I)=ROOTX
         !
130      CONTINUE
         TRIAL=(ROOTL+EV(I))*HALF
         !##IF CRAY (oldcode_for_cray)
         !            IF(TRIAL.EQ.ROOTL.OR.TRIAL.EQ.EV(I)) GOTO 200
         !##ELSE (oldcode_for_cray)
         EVDIFF=ABS(ROOTL-EV(I))
         IF(EVDIFF.LT.THETA1) GOTO 200
         IF(EVDIFF*ETAR.LT.ABS(TRIAL)) GOTO 200
         !##ENDIF (oldcode_for_cray)
         NOMTCH=N
         J=1
150      CONTINUE
         F0=A(J)-TRIAL
160      CONTINUE
         IF(ABS(F0).LT.THETA1) GOTO 170
         IF(F0.GE.ZERO) NOMTCH=NOMTCH-1
         J=J+1
         IF(J.GT.N) GOTO 180
         F0=A(J)-TRIAL-B(J-1)/F0
         GOTO 160
170      CONTINUE
         J=J+2
         NOMTCH=NOMTCH-1
         IF(J.LE.N) GOTO 150
180      CONTINUE
         IF(NOMTCH.GE.I) GOTO 190
         ROOTL=TRIAL
         GOTO 130
190      CONTINUE
         EV(I)=TRIAL
         NOM=MIN(NEVADD,NOMTCH)
         !            NOM=MIN0(NEVADD,NOMTCH)   		!original 
         EV(NOM)=TRIAL
         GOTO 130
200      CONTINUE
      ENDDO
      !
      ! Finished computing requested eigenvalues
      DO I=1,NEV
         EV(I)=EV(I+NADD)
      ENDDO
      !
      ! Compute eigenvectors (backtransformation)
      DO I=1,NEV
         AROOT=EV(I)
         DO J=1,N
            Y(J)=ONE
         ENDDO
         IA=IA+1
         IF(I.EQ.1) THEN
            IA=0
         ELSE
            IF(ABS(EV(I-1)-AROOT).GE.TOLER) IA=0
         ENDIF
         ELIM1=A(1)-AROOT
         ELIM2=W(1)
         DO J=1,NN
            IF(ABS(ELIM1).LE.ABS(W(J))) THEN
               TA(J)=W(J)
               TB(J)=A(J+1)-AROOT
               P(J)=W(J+1)
               TEMP=ONE
               IF(ABS(W(J)).GT.THETA1) TEMP=ELIM1/W(J)
               ELIM1=ELIM2-TEMP*TB(J)
               ELIM2=-TEMP*W(J+1)
            ELSE
               TA(J)=ELIM1
               TB(J)=ELIM2
               P(J)=ZERO
               TEMP=W(J)/ELIM1
               ELIM1=A(J+1)-AROOT-TEMP*ELIM2
               ELIM2=W(J+1)
            ENDIF
            B(J)=TEMP
         ENDDO
         !
         TA(N)=ELIM1
         TB(N)=ZERO
         P(N)=ZERO
         P(NN)=ZERO
         ITER=1
         IF(IA.NE.0) GOTO 460
         !
320      L=N+1
         DO J=1,N
            L=L-1
330         CONTINUE
            IF(L.EQ.N) THEN
               ELIM1=Y(L)
            ELSE IF(L.EQ.N-1) THEN
               ELIM1=Y(L)-Y(L+1)*TB(L)
            ELSE
               ELIM1=Y(L)-Y(L+1)*TB(L)-Y(L+2)*P(L)
            ENDIF
            !
            ! Overflow check
            IF(ABS(ELIM1).GT.DELBIG) THEN
               DO K=1,N
                  Y(K)=Y(K)/DELBIG
               ENDDO
               GOTO 330
            ENDIF
            TEMP=TA(L)
            IF(ABS(TEMP).LT.DELTA) TEMP=DELTA
            Y(L)=ELIM1/TEMP
         ENDDO
         !
         IF(ITER.EQ.2) GOTO 500
         ITER=ITER+1
         !
420      CONTINUE
         ELIM1=Y(1)
         DO J=1,NN
            IF(TA(J).EQ.W(J)) THEN
               Y(J)=Y(J+1)
               ELIM1=ELIM1-Y(J+1)*B(J)
            ELSE
               Y(J)=ELIM1
               ELIM1=Y(J+1)-ELIM1*B(J)
            ENDIF
         ENDDO
         Y(N)=ELIM1
         GOTO 320
         !
460      CONTINUE
         DO J=1,N
            RAND1=MOD(4099.0*RAND1,RPOWER)
            Y(J)=RAND1/RPOW1-ONE
         ENDDO
         GOTO 320
         !
         ! Orthog to previous
500      IF(IA.EQ.0) GOTO 550
         DO J1=1,IA
            K=I-J1
            TEMP=ZERO
            IPT=(K-1)*N
            DO J=1,N
               IPT=IPT+1
               TEMP=TEMP+Y(J)*VEC(IPT)
            ENDDO
            IPT=(K-1)*N
            DO J=1,N
               IPT=IPT+1
               Y(J)=Y(J)-TEMP*VEC(IPT)
            ENDDO
         ENDDO
550      CONTINUE
         IF(ITER.EQ.1) GOTO 420
         !
         ! Normalize
560      CONTINUE
         ELIM1=ZERO
         DO J=1,N
            ELIM1=MAX(ELIM1,ABS(Y(J)))
         ENDDO
         TEMP=ZERO
         DO J=1,N
            ELIM2=Y(J)/ELIM1
            TEMP=TEMP+ELIM2*ELIM2
         ENDDO
         TEMP=ONE/(SQRT(TEMP)*ELIM1)
         DO J=1,N
            Y(J)=Y(J)*TEMP
            IF(ABS(Y(J)).LT.DEL1) Y(J)=ZERO
         ENDDO
         IPT=(I-1)*N
         DO J=1,N
            IPT=IPT+1
            VEC(IPT)=Y(J)
         ENDDO
      ENDDO
      !   
      DO I=1,NEV
         IPT=(I-1)*N
         DO J=1,N
            IPT=IPT+1
            Y(J)=VEC(IPT)
         ENDDO
         !
         L=N-2
         MK=(N*(N-1))/2-3
         MK1=3
         !
         DO J=1,L
            T=ZERO
            K=N-J-1
            M=K+1
            DO KK=M,N
               T=T+DD(MK+KK)*Y(KK)
            ENDDO
            DO KK=M,N
               EPR=T*DD(MK+KK)
               Y(KK)=Y(KK)-EPR-EPR
            ENDDO
            MK=MK-MK1
            MK1=MK1+1
         ENDDO
         !
         T=ZERO
         DO J=1,N
            T=T+Y(J)*Y(J)
         ENDDO
         XNORM=SQRT(T)
         XNORM1=ONE/XNORM
         DO J=1,N
            Y(J)=Y(J)*XNORM1
         ENDDO
         !
         IPT=(I-1)*N
         DO J=1,N
            IPT=IPT+1
            VEC(IPT)=Y(J)
         ENDDO
      ENDDO
      !
      DO I=1,N
         EV(I)=EV(I)*ANORM
      ENDDO
      !
      RETURN
    END SUBROUTINE CHARMM_DIAGQ

  subroutine err_stop(info,id)
    !
    ! PROGRAM STOPS AT ERROR
    !
    implicit none
    character*(*),intent(in)  :: info
    integer,optional,intent(in)   :: id
    
    if(.not. present(id)) write(*,'(A,A)') 'ERRSTP> '//info
    if(present(id).and.id<1E08) write(*,'(A,I8)') 'ERRSTP> '//info,id
    if(present(id).and.id>=1E08) write(*,'(A,I12)') 'ERRSTP> '//info,id
    stop
    
    return
  end subroutine err_stop
    !
end module dpa
Program Multipe_chain_Dynamics_Perturbation_Analysis
!c
!c
!c       wcc: scaling interactions between CA -- CA
!c       wlg: scaling interactions between CA -- Ligand
!c       wll: scaling interactions between Ligand -- Ligand
!c       intwcc: gamma for interactions between two chains
!c       cutcc: cutoff for interactions between CA -- CA
!c       cutlg: cutoff for interactions between CA -- Ligand
!c       cutll: cutoff for interactions between Ligand -- Ligand
!c       intcutcc: cutoff for interactions between two chains
!        !nmodes: number of modes used for DPA calc.
!c
  implicit none
  character*100  file_cacrd,file_surfp,file_dpa,infocalc
  integer        calc_info	
  integer        nca,nsurf
  real*8         wcc,wctc,wlg,wll,cutcc,cutlg,cutll
  real*8         intwcc,intcutcc
  integer        nmodes
  real*8         modekappa
  real*8         ev_threshold	!maximum positive value for zero-mode eigenvalue
  real   time1,time2
!c
  read(*,*) file_cacrd,file_surfp,file_dpa
  read(*,*) wcc,wctc,wlg,wll,intwcc
  read(*,*) cutcc,cutlg,cutll,intcutcc
  read(*,*) nmodes,modekappa                !control diagonalization size
  read(*,*) ev_threshold
!c       
  call initial(file_cacrd,nca) 
  call initial(file_surfp,nsurf)
!c
  if(nmodes.le.0.or.nmodes.gt.nca*3) nmodes=nca*3
  call gene_surfdpa(calc_info,infocalc,ev_threshold,&
       file_cacrd,file_surfp,file_dpa,nca,nsurf,&
       cutcc,cutlg,cutll,wcc,wctc,wlg,wll,intcutcc,intwcc,&
       nmodes,modekappa,&
       time1,time2)

  if(calc_info.eq.0) then
     write(*,'(1x,I5,3x,F10.3)') calc_info,time2-time1
  else
     write(*,'(1x,I5,3x,A)') calc_info, infocalc
  endif
!
  stop
end Program Multipe_chain_Dynamics_Perturbation_Analysis

subroutine gene_surfdpa(calc_info,infocalc,ev_threshold,&
            file_cacrd,file_surfp,file_dpa,nca,nsurf,&
            cutcc,cutlg,cutll,wcc,wctc,wlg,wll,intcutcc,intwcc,&
            nmodes,modekappa,&
            time1,time2)
  implicit none
  character*100 infocalc,file_cacrd,file_surfp,file_dpa
  integer   calc_info,nca,nsurf,nca3,ndimca,ndimhet
  INTEGER   nat0,nat30,nhet,nhet3,nmodes
  real*8    modekappa
  PARAMETER (nat0=5000,nat30=3*nat0,nhet=1,nhet3=3*nhet)
  real*8    cacrd(3*nca),surfcrd(3*nsurf),hetcrd(nhet3)
  character*1  atmchain(nca),tmpresname*3
  real*8    wcc,wctc,wlg,wll,cutcc,cutlg,cutll,intcutcc,intwcc,pert
  real*8    ev_threshold
  real      time1,time2
  integer   i,k
!c
  nca3=3*nca
  ndimca=(nca3*(nca3+1))/2
  ndimhet=(nhet3*(nhet3+1))/2
!c       
  call SECOND(time1)
  open(unit=20,file=file_cacrd,STATUS='OLD',Err=5000) !Protein Backbone Coordinates
  do i=1,nca
     read(20,*,Err=5010) cacrd(3*i-2),cacrd(3*i-1),cacrd(3*i),tmpresname,atmchain(i)
  enddo
  close(20)
  open(unit=20,file=file_surfp,STATUS='OLD',Err=5500) !Protein Backbone Coordinates
  do i=1,nsurf
     read(20,*,Err=5510) surfcrd(3*i-2),surfcrd(3*i-1),surfcrd(3*i)
  enddo
  close(20)
!c
  call hessian_apoprotein(calc_info,infocalc,cacrd,atmchain,nca,nca3,ndimca,&
       cutcc,wcc,wctc,intcutcc,intwcc,ev_threshold,nmodes)
  if(calc_info.ne.0) return
  !c
  open(20,file=file_dpa,status='unknown')
  do i=1,nsurf
     hetcrd(1)=surfcrd(3*i-2)
     hetcrd(2)=surfcrd(3*i-1)
     hetcrd(3)=surfcrd(3*i)
     call DPA_1th_perturb(calc_info,infocalc,cacrd,hetcrd,nca+nhet,nca3,nhet3,&
          ndimca,ndimhet,cutlg,cutll,wlg,wll,nmodes,ev_threshold,pert)
     if(calc_info.ne.0) return
     write(20,'(1x,3F9.3,2x,F12.6,I10)') (hetcrd(k),k=1,3),pert,i
  enddo
  close(20)
!c
  call SECOND(time2)
!c
  calc_info=0
  return
!c
5000 calc_info=-9999
  infocalc='open cacrd file Error'
  return
5010 calc_info=-9998
  infocalc='read cacrd Error'
  return
5500 calc_info=-9997
  infocalc='open surfcrd file Error'
  return
5510 calc_info=-9996
  infocalc='read surfcrd Error'
  return
!
end subroutine gene_surfdpa

subroutine initial(file_input,n_point)
!c       
!c       return the number of point in file: file_input
!c       
  implicit none
  integer   n_point
  character*100 file_input
!c       
  open(unit=20,file=file_input,status='old',err=100)
  n_point=0
10 read(20,*,err=100,end=200) 
  n_point=n_point+1
  goto 10
50 write(*,*) -1, '  error occur in OPEN input crd file:',file_input
  stop 
100 write(*,*) -1, '  error occur in read input crd file:',file_input
  stop 
200 close(20)
!c       
  return
end subroutine initial

subroutine hessian_apoprotein(calc_info,infocalc,cacrd,atmchain,nca,nca3,ndimca,&
     cutcc,wcc,wctc,intcutcc,intwcc,ev_threshold,nmodes)
  use dpa
  implicit none
  character infocalc*100
  integer   calc_info,nca,nca3,ndimca
  real*8    cacrd(*),x(nca),y(nca),z(nca)
  character*1  atmchain(*)
  real*8 hessian1,hessian2,hessian3
  real*8 hessian4,hessian5,hessian6
  real*8 xx,yy,zz,sij,dsqrt,rr2ij
  real*8 cutoff,cutcc,wcc,wctc,intcutcc,intwcc,wiwj
  integer nmodes
  real*8  ev_threshold
  real*8,dimension(:),allocatable  :: hessian
  integer i,j,k,i1,j1,ki0,kj0,kh,isallocated
  !
  !
  allocate(hessian(ndimca),stat=isallocated)          !reate atom link for a residue
  if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR APOPROT HESSIAN(N)');
  hessian=0.0
!c   
!c--     Reading CA 
  do i=1,nca
     x(i)=cacrd(3*i-2)
     y(i)=cacrd(3*i-1)
     z(i)=cacrd(3*i)
  enddo
!c       
!c--  GENERATE HESSIAN
!c	!Off-diagonalize elements for lower-triangle off-diagonal blocks.
!c       i1=4 !(4,1) is the 1st element of the 2nd block of 1st column. 
!c       j1=1
!c       |1 2 3|
!c       |  4 5| main diagonal block
!c       |    6|
!c       
!c
  do i=1,nca		!atom i
     do j=i+1,nca	                 !atom j
        i1=3*(i-1)+1
        ki0=(i1-1)*nca3-((i1-1)*(i1-2))/2
        j1=3*(j-1)+1
        xx=x(j)-x(i)
        yy=y(j)-y(i)
        zz=z(j)-z(i)
        sij=sqrt(xx*xx+yy*yy+zz*zz)    !distance between  atom I and atom J
        rr2ij=1.0/sij/sij
        if(sij.lt.0.3) then
           calc_info=-5000
           infocalc='  Input Protein CRD error: Pair distance < 0.3'
           return
        endif
        !determin CUTOFF & INTERACTION strengh
        if(abs(abs(i-j)-1.).le.1.D-5) then
           wiwj=wctc
           ! write(*,*) i,j,wctc,wcc
        else
           wiwj=wcc
        endif
        cutoff=cutcc
        if(atmchain(i).ne.atmchain(j)) then
           wiwj=intwcc
           cutoff=intcutcc
        endif
!c       
        if(sij.le.cutoff) then
!c
           hessian1=-xx*xx*rr2ij*wiwj
           hessian2=-yy*xx*rr2ij*wiwj
           hessian3=-zz*xx*rr2ij*wiwj
           hessian4=-yy*yy*rr2ij*wiwj
           hessian5=-zz*yy*rr2ij*wiwj
           hessian6=-zz*zz*rr2ij*wiwj
           i1=3*(i-1)+1
           j1=3*(j-1)+1
           ki0=(i1-1)*NCA3-((i1-1)*(i1-2))/2
           kh=ki0+j1-i1+1 !kh for "line i1, column j1"
           hessian(kh)=hessian1 !line i1, column j1
           hessian(kh+1)=hessian2 !line i1, column j1+1
           hessian(kh+2)=hessian3 !line i1, column j1+2
!c       
           kh=ki0+NCA3-i1+1+j1-i1 !kh for "line i1+1, column j1"
           hessian(kh)=hessian2 !line i1+1, column j1
           hessian(kh+1)=hessian4 !line i1+1, column j1+1
           hessian(kh+2)=hessian5 !line i1+1, column j1
!c       
           kh=ki0+NCA3-i1+1+NCA3-i1+j1-i1-1 !kh for "line i1+2, column j1"
           hessian(kh)=hessian3 !line i1+2, column j1
           hessian(kh+1)=hessian5 !line i1+2, column j1
           hessian(kh+2)=hessian6 !line i1+2, column j1   
        endif
     enddo
  enddo

! --Diagonal-Block of Hessian
  do i=1,nca
     i1=(i-1)*3+1 
     ki0=(i1-1)*NCA3-((i1-1)*(i1-2))/2
     hessian1=0.0
     hessian2=0.0
     hessian3=0.0
     hessian4=0.0
     hessian5=0.0
     hessian6=0.0
     do j=i+1,nca
        j1=3*(j-1)+1
        kh=ki0+j1-i1+1	!kh for "line i1, column j1"
        hessian1=hessian1+hessian(kh) !line i1, column j1
        hessian2=hessian2+hessian(kh+1) !line i1, column j1+1
        hessian3=hessian3+hessian(kh+2) !line i1, column j1+2
        
        kh=ki0+NCA3-i1+1+j1-i1 !kh for "line i1+1, column j1"
        hessian4=hessian4+hessian(kh+1) !line i1+1, column j1+1
        hessian5=hessian5+hessian(kh+2) !line i1+1, column j1+2
!c       
        kh=ki0+NCA3-i1+1+NCA3-i1+j1-i1-1 !kh for "line i1+2, column j1"
        hessian6=hessian6+hessian(kh+2) !line i1+2, column j1+2
     enddo
     do j=1,i-1		!set j as row index now.
        j1=3*(j-1)+1
        kj0=(j1-1)*NCA3-((j1-1)*(j1-2))/2
!c       
        kh=kj0+i1-j1+1	!kh for "line j1, column i1"
        hessian1=hessian1+hessian(kh) !line j1, column i1
        hessian2=hessian2+hessian(kh+1) !line j1, column i1+1
        hessian3=hessian3+hessian(kh+2) !line j1, column i1+2
!c       
        kh=kj0+NCA3-j1+1+i1-j1 !kh for "line j1+1, column i1"
        hessian4=hessian4+hessian(kh+1) !line j1, column i1+1
        hessian5=hessian5+hessian(kh+2) !line j1, column i1+2
!c       
        kh=kj0+NCA3-j1+1+NCA3-j1+i1-j1-1 !kh for "line j1+2, column i1"
        hessian6=hessian6+hessian(kh+2) !line j1, column i1+2  
!c     
     enddo
     kh=ki0+1
     hessian(kh)=-hessian1
     hessian(kh+1)=-hessian2
     hessian(kh+2)=-hessian3
     kh=ki0+NCA3-i1+1
     hessian(kh+1)=-hessian4
     hessian(kh+2)=-hessian5
     kh=ki0+NCA3-i1+1+NCA3-i1
     hessian(kh+1)=-hessian6    
  enddo
!c
!  call CharmmDiagq(nca3,ione,nmodes,ev0,vec0)
  call cdiag(nca3,nmodes,hessian,ev0,vec0)
  if(dabs(ev0(6)).gt.ev_threshold.or.ev0(7).lt.ev_threshold) then 
     calc_info=-1
     infocalc='  CUTCC for protein is too small'
     return
  endif
!c
  calc_info=0
!
  return
end subroutine hessian_apoprotein

subroutine DPA_1th_perturb(calc_info,infocalc,cacrd,hetcrd,&
  nsys,nca3,nhet3,ndimca,ndimhet,cutlg,cutll,wlg,wll,nmodes,&
  ev_threshold,lgpert)
  use dpa
  implicit none
  character*100 infocalc
  integer   calc_info,nsys,nca,nhet,nca3,nhet3,ndimca,ndimhet,nmodes
  real*8    cacrd(*),hetcrd(*)
  real*8    cutoff,cutlg,cutll,wlg,wll,wiwj
!
  integer   nat1,nat31 !nat1: maximum number of atoms that ligand binds
  parameter (nat1=200,nat31=3*nat1)

  real*8    ev_threshold
  real*8    gessian(nca3,nhet3),tmpw(nhet3),g_sqrt_invk(nca3,nhet3)
  real*8,dimension(:),allocatable :: kessian,kev,kvec
  real*8    ev(nca3)
  integer   rp(nca3),idgnz(nca3),col_id(nca3*nat31),h_id(nca3*nat31)
  real*8    pessian(nca3*nat31),hv(nca3)
  real*8    x(nsys),y(nsys),z(nsys),xx,yy,zz,sij,dsqrt,rr2ij
  real*8    hessian1,hessian2,hessian3,hessian4,hessian5,hessian6
  real*8    lgpert,tmp
  integer   i,j,k,l,i0,j0,i1,j1,k1,ki0,kj0,kh,isallocated
!c
  nca=nca3/3
  nhet=nhet3/3
  allocate(kessian(ndimhet),stat=isallocated)          !reate atom link for a residue
  if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR KESSIAN(N)');
  kessian=0.0
!c   
!c--     Reading CA & HETATM coordinates
  do i=1,nca
     x(i)=cacrd(3*i-2)
     y(i)=cacrd(3*i-1)
     z(i)=cacrd(3*i)
  enddo
  do i=1,nhet
     x(nca+i)=hetcrd(3*i-2)
     y(nca+i)=hetcrd(3*i-1)
     z(nca+i)=hetcrd(3*i)
  enddo
!c       
!c--  GENERATE HESSIAN
!c	!Off-diagonalize elements for lower-triangle off-diagonal blocks.
!c       i1=4 !(4,1) is the 1st element of the 2nd block of 1st column. 
!c       j1=1
!c       |1 2 3|
!c       |  4 5| main diagonal block
!c       |    6|
!c       
  do i=1,nca3
     idgnz(i)=0		!idex for NONZERO elements in GESSIAN
  enddo
  do i=1,nca3
     do j=1,nhet3
        gessian(i,j)=0.0
     enddo
  enddo
  do i=1,ndimhet
     kessian(i)=0.0
  enddo
!c
  do i=1,nsys-1		!i<=nca: Protein residues; i>nca: Ligand residues
     do j=nca+1,nsys	!ligand residues
        xx=x(j)-x(i)
        yy=y(j)-y(i)
        zz=z(j)-z(i)
        sij=sqrt(xx*xx+yy*yy+zz*zz) !distance between  atom I and atom J
        rr2ij=1.0/sij/sij
        if(sij.lt.1.) then
           calc_info=-5100
           infocalc='  Input CRD data Error: Protein/Ligand distance < 1' 
           return
        endif
        !determin CUTOFF & INTERACTION strengh
        if(i.le.nca) then !i: protein/ligand interaction :: protein atom; j: ligand
           wiwj= wlg
           cutoff=cutlg
        else		!i,j: interaction within ligand atoms
           wiwj=wll
           cutoff=cutll
        endif
!c            
        if(sij.le.cutoff) then
           hessian1=-xx*xx*rr2ij*wiwj
           hessian2=-yy*xx*rr2ij*wiwj
           hessian3=-zz*xx*rr2ij*wiwj
           hessian4=-yy*yy*rr2ij*wiwj
           hessian5=-zz*yy*rr2ij*wiwj
           hessian6=-zz*zz*rr2ij*wiwj
           if (i.le.nca) then ! For GESSIAN matrix -- Protein/Ligand Interaction portion
              i1=3*(i-1)+1
              j1=3*(j-nca-1)+1
              gessian(i1,j1)=hessian1     !line i1,   column j1
              gessian(i1,j1+1)=hessian2   !line i1,   column j1+1
              gessian(i1,j1+2)=hessian3   !line i1,   column j1+2
              gessian(i1+1,j1)=hessian2   !line i1+1, column j1
              gessian(i1+1,j1+1)=hessian4 !line i1+1, column j1+1
              gessian(i1+1,j1+2)=hessian5 !line i1+1, column j1+2
              gessian(i1+2,j1)=hessian3   !line i1+2, column j1
              gessian(i1+2,j1+1)=hessian5 !line i1+2, column j1+1
              gessian(i1+2,j1+2)=hessian6 !line i1+2, column j1+2  
              idgnz(i1)=1
              idgnz(i1+1)=1
              idgnz(i1+2)=1
           else		! For KESSIAN matrix -- Ligand portion
              i1=3*(i-nca-1)+1
              ki0=(i1-1)*nhet3-((i1-1)*(i1-2))/2
              j1=3*(j-nca-1)+1		    
              
              kh=ki0+j1-i1+1 !kh for "line i1, column j1"
              kessian(kh)=hessian1 !line i1, column j1
              kessian(kh+1)=hessian2 !line i1, column j1+1
              kessian(kh+2)=hessian3 !line i1, column j1+2
!c       
              kh=ki0+nhet3-i1+1+j1-i1 !kh for "line i1+1, column j1"
              kessian(kh)=hessian2 !line i1+1, column j1
              kessian(kh+1)=hessian4 !line i1+1, column j1+1
              kessian(kh+2)=hessian5 !line i1+1, column j1
!c       
              kh=ki0+nhet3-i1+1+nhet3-i1+j1-i1-1 !kh for "line i1+2, column j1"
              kessian(kh)=hessian3 !line i1+2, column j1
              kessian(kh+1)=hessian5 !line i1+2, column j1
              kessian(kh+2)=hessian6 !line i1+2, column j1   
           endif
        endif
     enddo
  enddo
!c
!diagonal-Block of Kessian
  do i=1,nhet
     i1=(i-1)*3+1 
     ki0=(i1-1)*nhet3-((i1-1)*(i1-2))/2
     hessian1=0.0
     hessian2=0.0
     hessian3=0.0
     hessian4=0.0
     hessian5=0.0
     hessian6=0.0
     do j=i+1,nhet
        j1=3*(j-1)+1
        kh=ki0+j1-i1+1	!kh for "line i1, column j1"
        hessian1=hessian1+Kessian(kh) !line i1, column j1
        hessian2=hessian2+Kessian(kh+1) !line i1, column j1+1
        hessian3=hessian3+Kessian(kh+2) !line i1, column j1+2
!c
        kh=ki0+nhet3-i1+1+j1-i1 !kh for "line i1+1, column j1"
        hessian4=hessian4+Kessian(kh+1) !line i1+1, column j1+1
        hessian5=hessian5+Kessian(kh+2) !line i1+1, column j1+2
!c       
        kh=ki0+nhet3-i1+1+nhet3-i1+j1-i1-1 !kh for "line i1+2, column j1"
        hessian6=hessian6+Kessian(kh+2) !line i1+2, column j1+2
     enddo
     do j=1,i-1		!set j as row index now.
        j1=3*(j-1)+1
        kj0=(j1-1)*nhet3-((j1-1)*(j1-2))/2
!c       
        kh=kj0+i1-j1+1	!kh for "line j1, column i1"
        hessian1=hessian1+Kessian(kh) !line j1, column i1
        hessian2=hessian2+Kessian(kh+1) !line j1, column i1+1
        hessian3=hessian3+Kessian(kh+2) !line j1, column i1+2
!c       
        kh=kj0+nhet3-j1+1+i1-j1 !kh for "line j1+1, column i1"
        hessian4=hessian4+Kessian(kh+1) !line j1, column i1+1
        hessian5=hessian5+Kessian(kh+2) !line j1, column i1+2
!c       
        kh=kj0+nhet3-j1+1+nhet3-j1+i1-j1-1 !kh for "line j1+2, column i1"
        hessian6=hessian6+Kessian(kh+2) !line j1, column i1+2  
     enddo
     do j=1,nca
        j1=3*(j-1)+1
        hessian1=hessian1+gessian(j1,i1)
        hessian2=hessian2+gessian(j1,i1+1)
        hessian3=hessian3+gessian(j1,i1+2)
        hessian4=hessian4+gessian(j1+1,i1+1)
        hessian5=hessian5+gessian(j1+1,i1+2)
        hessian6=hessian6+gessian(j1+2,i1+2)
     enddo
     kh=ki0+1
     kessian(kh)=-hessian1
     kessian(kh+1)=-hessian2
     kessian(kh+2)=-hessian3
     kh=ki0+nhet3-i1+1
     kessian(kh+1)=-hessian4
     kessian(kh+2)=-hessian5
     kh=ki0+nhet3-i1+1+nhet3-i1
     kessian(kh+1)=-hessian6
  enddo

  call cdiag(nhet3,nhet3,kessian,kev,kvec)

  do i=1,nhet3		!check inver of  Kessian        
     if(Dabs(kev(i)) .le. 1.d-8) then
        calc_info=-9
        infocalc=' : Failure: SINGULAR MATRIX FOR KESSIAN'
        return 
     endif
  enddo
!c
!c-- 2)generate quasi inverse of Kessian
  do i=1,nhet3
     tmpw(i)=1./dsqrt(kev(i))
  enddo
  j0=0
  do j=1,nhet3
     do i=1,nhet3
        kvec(j0+i)=kvec(j0+i)*tmpw(j)
     enddo
     j0=j0+nhet3
  enddo
!c
!c--     3) generate qausi-perturbed Matrix (QpM): Gessian x SQRT(Kessian^-1) = GESS x KESS^(-1/2)
  do i=1,nca
     do j=1,nhet3
        g_sqrt_invk(i,j)=0.
     enddo
  enddo
  do i=1,nca3
     if(idgnz(i).ne.0) then
        j0=0
        do j=1,nhet3
           tmp=0.
           do k=1,nhet3
              tmp=tmp+gessian(i,k)*kvec(j0+k)
           enddo
           g_sqrt_invk(i,j)=tmp
           j0=j0+nhet3
        enddo
     endif
  enddo
!c--
  do i=1,nca3
     rp(i)=0		!column index of NONZERO element in ith row of  Perturbation Matrix 
  enddo
!c
  kh=0
  i1=0
  do i=1,nca3		!pessian(i,i)=0 if ith atom do NOT inter with Ligands
     j1=i1
     do j=i,nca3
        if(idgnz(j).ne.0.and.idgnz(i).ne.0) then
!c
!c-- CALC. PESSIAN
           tmp=0.
           do k=1,nhet3
              tmp=tmp+g_sqrt_invk(i,k)*g_sqrt_invk(j,k) !from  "G x (K^-1) x G"
           enddo
           if(j.le.(i+mod(3-mod(i,3),3))) then !Diagonal Term as required from sum_i_k=0
              i0=mod(j,3)
              if(i0.eq.0) i0=3
              do j0=i0,nhet3,3 
                 tmp=tmp+gessian(i,j0)
              enddo
           endif
           kh=kh+1
           rp(i)=rp(i)+1
           if(rp(i).gt.nat31) then
              calc_info=-2000
              infocalc=' Declared Dimension of RP (neighboring contacts) is too small'
              return
           endif
           pessian(kh)=-tmp 
!c
!c-- INDEX for PESSIAN
!c		 col_id(i,rp(i))=j !in row i, column index of NONZERO element is j in origianl Hessian
!c		 h_id(i,rp(i))=kh !in row i, nonzero element's value is the kh element of New Hessian
           col_id(i1+rp(i))=j
           h_id(i1+rp(i))=kh
           !c
           if(j.ne.i)  then
              rp(j)=rp(j)+1
              if(rp(j).gt.nat31) then
                 calc_info=-2000
                 infocalc=' Declared Dimension of RP (neighboring contacts) is too small'
                 return
              endif
!c		    h_id(j,rp(j))=kh !in row j, nonzero element's value is kh_th element of New Hessian
!c		    col_id(j,rp(j))=i !in row j,column index of NONZERO element is i_th in origianl Hessian
              col_id(j1+rp(j))=i
              h_id(j1+rp(j))=kh
!c
           endif
        endif
        j1=j1+nat31
     enddo
     i1=i1+nat31
  enddo
!c
  i0=0
!c	do i=1,nca3
  do i=1,nmodes  !instead of using all 3N modes, using fewer low-frequency modes
     k1=0
     do k=1,nca3
        tmp=0.
        do j=1,rp(k)
!c		 tmp=tmp+pessian(h_id(k,j))*vec0(i0+col_id(k,j)) !Sum_j {P(k,j)*V(j,i)}
           tmp=tmp+pessian(h_id(k1+j))*vec0(i0+col_id(k1+j)) !Sum_j {P(k,j)*V(j,i)}
        enddo
        hv(k)=tmp		!H(k,i)
        k1=k1+nat31
     enddo
     tmp=0.
     do k=1,nca3
        tmp=tmp+vec0(k+i0)*hv(k) !Sum_k {V(k,i)*H(k,i)}
     enddo
     ev(i)=ev0(i)+tmp	!write(11,'(1x,I8,2x,3F12.5)') i,tmp,ev(i),ev0(i)
     i0=i0+nca3
  enddo
!c	
  lgpert=0.0
!c	do i=7,nca3
  do i=7,nmodes   !instead of using all 3N modes, using fewer low-frequency modes
     if(ev(i).lt.ev_threshold) then
        lgpert=-9999.
        return
     else
        lgpert=lgpert+dlog(ev(i)/ev0(i))
     endif
  enddo
  calc_info=0
!c
  lgpert=0.0
!c	do i=7,nca3
  do i=7,nmodes   !instead of using all 3N modes, using fewer low-frequency modes
     if(ev(i).lt.ev_threshold) then
        lgpert=-9999.
        return
     else
        lgpert=lgpert+dlog(ev(i)/ev0(i))
     endif
  enddo
  calc_info=0
!c
  return
end subroutine DPA_1th_perturb
