source/unres/src-HCD-5D/econstrq.F

   1       subroutine EconstrQ
   2 c     MD with umbrella_sampling using Wolyne's distance measure as a constraint
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.VAR'
   7       include 'COMMON.MD'
   8 #ifndef LANG0
   9       include 'COMMON.LANGEVIN'
  10 #else
  11       include 'COMMON.LANGEVIN.lang0'
  12 #endif
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.DERIV'
  15       include 'COMMON.GEO'
  16       include 'COMMON.LOCAL'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.NAMES'
  20       include 'COMMON.TIME1'
  21       double precision uzap1,uzap2,hm1,hm2,hmnum
  22       double precision ucdelan,dUcartan(3,0:MAXRES)
  23      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
  24      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
  25       integer kstart,kend,lstart,lend,idummy
  26       double precision delta /1.0d-7/
  27       do i=0,nres
  28          do j=1,3
  29             duconst(j,i)=0.0d0
  30             dudconst(j,i)=0.0d0
  31             duxconst(j,i)=0.0d0
  32             dudxconst(j,i)=0.0d0
  33          enddo
  34       enddo
  35       Uconst=0.0d0
  36 #ifdef DEBUG
  37       write (iout,*) "econstrq: nfrag",nfrag," iset",iset,
  38      &   " ifrag",(ifrag(1,i,iset),ifrag(2,i,iset),i=1,nfrag)
  39       do i=1,nfrag
  40          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
  41      &    ,idummy,idummy)
  42          write (iout,*) "fragment",i," qfrag",qfrag(i)
  43 #endif
  44          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
  45 c Calculating the derivatives of Constraint energy with respect to Q
  46          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
  47      &     qinfrag(i,iset))
  48 c         hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
  49 c        hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
  50 c         hmnum=(hm2-hm1)/delta
  51 c         write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
  52 c     &   qinfrag(i,iset))
  53 c         write(iout,*) "harmonicnum frag", hmnum
  54 c Calculating the derivatives of Q with respect to cartesian coordinates
  55          call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
  56      &   ,idummy,idummy)
  57 c         write(iout,*) "dqwol "
  58 c         do ii=1,nres
  59 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
  60 c         enddo
  61 c         write(iout,*) "dxqwol "
  62 c         do ii=1,nres
  63 c           write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
  64 c         enddo
  65 c Calculating numerical gradients of dU/dQi and dQi/dxi
  66 c        call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
  67 c     &  ,idummy,idummy)
  68 c  The gradients of Uconst in Cs
  69          do ii=0,nres
  70             do j=1,3
  71                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
  72                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
  73             enddo
  74          enddo
  75       enddo
  76       do i=1,npair
  77          kstart=ifrag(1,ipair(1,i,iset),iset)
  78          kend=ifrag(2,ipair(1,i,iset),iset)
  79          lstart=ifrag(1,ipair(2,i,iset),iset)
  80          lend=ifrag(2,ipair(2,i,iset),iset)
  81          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
  82          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
  83 c  Calculating dU/dQ
  84          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
  85 c         hm1=harmonic(qpair(i),qinpair(i,iset))
  86 c        hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
  87 c         hmnum=(hm2-hm1)/delta
  88 c         write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
  89 c     &   qinpair(i,iset))
  90 c         write(iout,*) "harmonicnum pair ", hmnum
  91 c Calculating dQ/dXi
  92          call qwolynes_prim(kstart,kend,.false.
  93      &   ,lstart,lend)
  94 c         write(iout,*) "dqwol "
  95 c         do ii=1,nres
  96 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
  97 c         enddo
  98 c         write(iout,*) "dxqwol "
  99 c         do ii=1,nres
 100 c          write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
 101 c        enddo
 102 c Calculating numerical gradients
 103 c        call qwol_num(kstart,kend,.false.
 104 c     &  ,lstart,lend)
 105 c The gradients of Uconst in Cs
 106          do ii=0,nres
 107             do j=1,3
 108                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 109                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 110             enddo
 111          enddo
 112       enddo
 113 c      write(iout,*) "Uconst inside subroutine ", Uconst
 114 c Transforming the gradients from Cs to dCs for the backbone
 115       do i=0,nres
 116          do j=i+1,nres
 117            do k=1,3
 118              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 119            enddo
 120          enddo
 121       enddo
 122 c  Transforming the gradients from Cs to dCs for the side chains
 123       do i=1,nres
 124          do j=1,3
 125            dudxconst(j,i)=duxconst(j,i)
 126          enddo
 127       enddo
 128 c      write(iout,*) "dU/ddc backbone "
 129 c       do ii=0,nres
 130 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 131 c      enddo
 132 c      write(iout,*) "dU/ddX side chain "
 133 c      do ii=1,nres
 134 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 135 c      enddo
 136 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 137 c      call dEconstrQ_num
 138       return
 139       end