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