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