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