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