source/unres/src-5hdiag-tmp/deconstrq_num.F

   1       subroutine dEconstrQ_num
   2 c Calculating numerical dUconst/ddc and dUconst/ddx
   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
  22       double precision dUcartan(3,0:MAXRES)
  23      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
  24       integer kstart,kend,lstart,lend,idummy
  25       double precision delta /1.0d-7/
  26 c     For the backbone
  27       do i=0,nres-1
  28          do j=1,3
  29             dUcartan(j,i)=0.0d0
  30             cdummy(j,i)=dc(j,i)
  31             dc(j,i)=dc(j,i)+delta
  32             call chainbuild_cart
  33             uzap2=0.0d0
  34             do ii=1,nfrag
  35              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
  36      &         ,idummy,idummy)
  37                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
  38      &          qinfrag(ii,iset))
  39             enddo
  40             do ii=1,npair
  41                kstart=ifrag(1,ipair(1,ii,iset),iset)
  42                kend=ifrag(2,ipair(1,ii,iset),iset)
  43                lstart=ifrag(1,ipair(2,ii,iset),iset)
  44                lend=ifrag(2,ipair(2,ii,iset),iset)
  45                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
  46                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
  47      &           qinpair(ii,iset))
  48             enddo
  49             dc(j,i)=cdummy(j,i)
  50             call chainbuild_cart
  51             uzap1=0.0d0
  52              do ii=1,nfrag
  53              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
  54      &         ,idummy,idummy)
  55                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
  56      &          qinfrag(ii,iset))
  57             enddo
  58             do ii=1,npair
  59                kstart=ifrag(1,ipair(1,ii,iset),iset)
  60                kend=ifrag(2,ipair(1,ii,iset),iset)
  61                lstart=ifrag(1,ipair(2,ii,iset),iset)
  62                lend=ifrag(2,ipair(2,ii,iset),iset)
  63                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
  64                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
  65      &          qinpair(ii,iset))
  66             enddo
  67             ducartan(j,i)=(uzap2-uzap1)/(delta)
  68          enddo
  69       enddo
  70 c Calculating numerical gradients for dU/ddx
  71       do i=0,nres-1
  72          duxcartan(j,i)=0.0d0
  73          do j=1,3
  74             cdummy(j,i)=dc(j,i+nres)
  75             dc(j,i+nres)=dc(j,i+nres)+delta
  76             call chainbuild_cart
  77             uzap2=0.0d0
  78             do ii=1,nfrag
  79              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
  80      &         ,idummy,idummy)
  81                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
  82      &          qinfrag(ii,iset))
  83             enddo
  84             do ii=1,npair
  85                kstart=ifrag(1,ipair(1,ii,iset),iset)
  86                kend=ifrag(2,ipair(1,ii,iset),iset)
  87                lstart=ifrag(1,ipair(2,ii,iset),iset)
  88                lend=ifrag(2,ipair(2,ii,iset),iset)
  89                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
  90                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
  91      &          qinpair(ii,iset))
  92             enddo
  93             dc(j,i+nres)=cdummy(j,i)
  94             call chainbuild_cart
  95             uzap1=0.0d0
  96              do ii=1,nfrag
  97                qfrag(ii)=qwolynes(ifrag(1,ii,iset),
  98      &          ifrag(2,ii,iset),.true.,idummy,idummy)
  99                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
 100      &          qinfrag(ii,iset))
 101             enddo
 102             do ii=1,npair
 103                kstart=ifrag(1,ipair(1,ii,iset),iset)
 104                kend=ifrag(2,ipair(1,ii,iset),iset)
 105                lstart=ifrag(1,ipair(2,ii,iset),iset)
 106                lend=ifrag(2,ipair(2,ii,iset),iset)
 107                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 108                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
 109      &          qinpair(ii,iset))
 110             enddo
 111             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 112          enddo
 113       enddo
 114       write(iout,*) "Numerical dUconst/ddc backbone "
 115       do ii=0,nres
 116         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 117       enddo
 118 c      write(iout,*) "Numerical dUconst/ddx side-chain "
 119 c      do ii=1,nres
 120 c         write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 121 c      enddo
 122       return
 123       end
 124 c---------------------------------------------------------------------------
 125