source/unres/src-HCD-5D/q_measure1.F

   1       double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.VAR'
   8       include 'COMMON.MD'
   9       integer i,j,jl,k,l,il,kl,nl,np,seg1,seg2,seg3,seg4,secseg
  10       integer nsep /3/
  11       double precision dist,qm
  12       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  13       logical lprn /.false./
  14       logical flag
  15       qq = 0.0d0
  16       nl=0
  17       do i=0,nres
  18         do j=1,3
  19           dqwol(j,i)=0.0d0
  20           dxqwol(j,i)=0.0d0
  21         enddo
  22       enddo
  23       if (lprn) then
  24       write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4,
  25      & " flag",flag
  26       call flush(iout)
  27       endif
  28       if (flag) then
  29         do il=seg1+nsep,seg2
  30           do jl=seg1,il-nsep
  31             nl=nl+1
  32             if (itype(il).ne.10) then
  33               ilnres=il+nres
  34             else
  35               ilnres=il
  36             endif
  37             if (itype(jl).ne.10) then
  38               jlnres=jl+nres
  39             else
  40               jlnres=jl
  41             endif
  42             qqijCM = qcontrib(il,jl,ilnres,jlnres)
  43             qq = qq+qqijCM
  44             if (lprn) then
  45               write (iout,*) "qqijCM",qqijCM
  46               call flush(iout)
  47             endif
  48           enddo
  49         enddo
  50         if (lprn) then
  51           write (iout,*) "nl",nl," qq",qq
  52           call flush(iout)
  53         endif
  54       else
  55         do il=seg1,seg2
  56           if((seg3-il).lt.3) then
  57              secseg=il+3
  58           else
  59              secseg=seg3
  60           endif
  61           do jl=secseg,seg4
  62             nl=nl+1
  63             if (itype(il).ne.10) then
  64               ilnres=il+nres
  65             else
  66               ilnres=il
  67             endif
  68             if (itype(jl).ne.10) then
  69               jlnres=jl+nres
  70             else
  71               jlnres=jl
  72             endif
  73             qqijCM = qcontrib(il,jl,ilnres,jlnres)
  74             qq = qq+qqijCM
  75             if (lprn) then
  76               write (iout,*) "qqijCM",qqijCM
  77               call flush(iout)
  78             endif
  79           enddo
  80         enddo
  81       endif
  82       qq = qq/nl
  83       qwolynes=1.0d0-qq
  84       do i=0,nres
  85         do j=1,3
  86           dqwol(j,i)=dqwol(j,i)/nl
  87           dxqwol(j,i)=dxqwol(j,i)/nl
  88         enddo
  89       enddo
  90       return
  91       end
  92 c-------------------------------------------------------------------
  93       subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
  94       implicit real*8 (a-h,o-z)
  95       include 'DIMENSIONS'
  96       include 'COMMON.IOUNITS'
  97       include 'COMMON.CHAIN'
  98       include 'COMMON.INTERACT'
  99       include 'COMMON.VAR'
 100       include 'COMMON.MD'
 101       integer seg1,seg2,seg3,seg4
 102       logical flag
 103       double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
 104      & qwolxan(3,0:maxres),q1,q2
 105       double precision delta /1.0d-7/
 106       write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4
 107       write(iout,*) "dQ/dc backbone "
 108        do i=0,nres
 109         write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3)
 110       enddo
 111       write(iout,*) "dQ/dX side chain "
 112       do i=1,nres
 113             write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3)
 114       enddo
 115       do i=1,nres
 116         do j=1,3
 117           cdummy(j,i)=c(j,i)
 118           c(j,i)=c(j,i)-delta
 119           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 120           c(j,i)=cdummy(j,i)+delta
 121           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 122           qwolan(j,i)=0.5d0*(q2-q1)/delta
 123           c(j,i)=cdummy(j,i)
 124 c          write (iout,*) "i",i," j",j," q1",q1," a2",q2
 125         enddo
 126       enddo
 127       do i=1,nres
 128         do j=1,3
 129           cdummy(j,i+nres)=c(j,i+nres)
 130           c(j,i+nres)=c(j,i+nres)-delta
 131           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 132           c(j,i+nres)=cdummy(j,i+nres)+delta
 133           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 134           qwolxan(j,i)=0.5d0*(q2-q1)/delta
 135           c(j,i+nres)=cdummy(j,i+nres)
 136         enddo
 137       enddo
 138       write(iout,*) "Numerical Q cartesian gradients backbone: "
 139       do i=0,nres
 140         write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
 141       enddo
 142       write(iout,*) "Numerical Q cartesian gradients side-chain: "
 143       do i=0,nres
 144         write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
 145       enddo
 146       return
 147       end
 148 c------------------------------------------------------------------------
 149       subroutine EconstrQ
 150 c     MD with umbrella_sampling using Wolyne's distance measure as a constraint
 151       implicit real*8 (a-h,o-z)
 152       include 'DIMENSIONS'
 153       include 'COMMON.CONTROL'
 154       include 'COMMON.VAR'
 155       include 'COMMON.MD'
 156 #ifndef LANG0
 157       include 'COMMON.LANGEVIN'
 158 #else
 159 #ifdef FIVEDIAG
 160       include 'COMMON.LANGEVIN.lang0.5diag'
 161 #else
 162       include 'COMMON.LANGEVIN.lang0'
 163 #endif
 164 #endif
 165       include 'COMMON.CHAIN'
 166       include 'COMMON.DERIV'
 167       include 'COMMON.GEO'
 168       include 'COMMON.LOCAL'
 169       include 'COMMON.INTERACT'
 170       include 'COMMON.IOUNITS'
 171       include 'COMMON.NAMES'
 172       include 'COMMON.TIME1'
 173       double precision uzap1,uzap2,hm1,hm2,hmnum
 174       double precision ucdelan,dUcartan(3,0:MAXRES)
 175      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
 176      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
 177       integer kstart,kend,lstart,lend,idummy
 178       double precision delta /1.0d-7/
 179       do i=0,nres
 180          do j=1,3
 181             duconst(j,i)=0.0d0
 182             dudconst(j,i)=0.0d0
 183             duxconst(j,i)=0.0d0
 184             dudxconst(j,i)=0.0d0
 185          enddo
 186       enddo
 187       Uconst=0.0d0
 188       do i=1,nfrag
 189          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 190      &   ,idummy,idummy)
 191          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
 192 c Calculating the derivatives of Constraint energy with respect to Q
 193          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
 194 c Calculating the derivatives of Q with respect to cartesian coordinates
 195          do ii=0,nres
 196             do j=1,3
 197                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
 198                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
 199             enddo
 200          enddo
 201 c      write (iout,*) "Calling qwol_num"
 202 c      call qwol_num(ifrag(1,i),ifrag(2,i),.true.,idummy,idummy)
 203       enddo
 204       do i=1,npair
 205          kstart=ifrag(1,ipair(1,i,iset),iset)
 206          kend=ifrag(2,ipair(1,i,iset),iset)
 207          lstart=ifrag(1,ipair(2,i,iset),iset)
 208          lend=ifrag(2,ipair(2,i,iset),iset)
 209          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
 210          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
 211 c  Calculating dU/dQ
 212          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
 213 c Calculating dQ/dXi
 214          do ii=0,nres
 215             do j=1,3
 216                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 217                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 218             enddo
 219          enddo
 220       enddo
 221 c      write(iout,*) "Uconst inside subroutine ", Uconst
 222 c Transforming the gradients from Cs to dCs for the backbone
 223       do i=0,nres
 224          do j=i+1,nres
 225            do k=1,3
 226              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 227            enddo
 228          enddo
 229       enddo
 230 c  Transforming the gradients from Cs to dCs for the side chains
 231       do i=1,nres
 232          do j=1,3
 233            dudxconst(j,i)=duxconst(j,i)
 234          enddo
 235       enddo
 236 c      write(iout,*) "dU/dc backbone "
 237 c       do ii=0,nres
 238 c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
 239 c      enddo
 240 c      write(iout,*) "dU/dX side chain "
 241 c      do ii=1,nres
 242 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 243 c      enddo
 244 c      write(iout,*) "dU/ddc backbone "
 245 c       do ii=0,nres
 246 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 247 c      enddo
 248 c      write(iout,*) "dU/ddX side chain "
 249 c      do ii=1,nres
 250 c            write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
 251 c      enddo
 252 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 253 c      call dEconstrQ_num
 254       return
 255       end
 256 c-----------------------------------------------------------------------
 257       subroutine dEconstrQ_num
 258 c Calculating numerical dUconst/ddc and dUconst/ddx
 259       implicit real*8 (a-h,o-z)
 260       include 'DIMENSIONS'
 261       include 'COMMON.CONTROL'
 262       include 'COMMON.VAR'
 263       include 'COMMON.MD'
 264 #ifndef LANG0
 265       include 'COMMON.LANGEVIN'
 266 #else
 267 #ifdef FIVEDIAG
 268       include 'COMMON.LANGEVIN.lang0.5diag'
 269 #else
 270       include 'COMMON.LANGEVIN.lang0'
 271 #endif
 272 #endif
 273       include 'COMMON.CHAIN'
 274       include 'COMMON.DERIV'
 275       include 'COMMON.GEO'
 276       include 'COMMON.LOCAL'
 277       include 'COMMON.INTERACT'
 278       include 'COMMON.IOUNITS'
 279       include 'COMMON.NAMES'
 280       include 'COMMON.TIME1'
 281       double precision uzap1,uzap2
 282       double precision dUcartan(3,0:MAXRES)
 283      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
 284       integer kstart,kend,lstart,lend,idummy
 285       double precision delta /1.0d-7/
 286 c     For the backbone
 287       do i=0,nres-1
 288          do j=1,3
 289             dUcartan(j,i)=0.0d0
 290             cdummy(j,i)=dc(j,i)
 291             dc(j,i)=dc(j,i)+delta
 292             call chainbuild_cart
 293             uzap2=0.0d0
 294             do ii=1,nfrag
 295                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 296      &         .true.,idummy,idummy)
 297                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
 298      &            qinfrag(ii,iset))
 299             enddo
 300             do ii=1,npair
 301                kstart=ifrag(1,ipair(1,ii,iset),iset)
 302                kend=ifrag(2,ipair(1,ii,iset),iset)
 303                lstart=ifrag(1,ipair(2,ii,iset),iset)
 304                lend=ifrag(2,ipair(2,ii,iset),iset)
 305                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 306                uzap2=uzap2+wpair(ii,iset)*
 307      &             harmonic(qpair(ii),qinpair(ii,iset))
 308             enddo
 309             dc(j,i)=cdummy(j,i)
 310             call chainbuild_cart
 311             uzap1=0.0d0
 312              do ii=1,nfrag
 313                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 314      &         .true.,idummy,idummy)
 315                uzap1=uzap1+wfrag(ii,iset)*
 316      &                 harmonic(qfrag(ii),qinfrag(ii,iset))
 317             enddo
 318             do ii=1,npair
 319                kstart=ifrag(1,ipair(1,ii,iset),iset)
 320                kend=ifrag(2,ipair(1,ii,iset),iset)
 321                lstart=ifrag(1,ipair(2,ii,iset),iset)
 322                lend=ifrag(2,ipair(2,ii,iset),iset)
 323                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 324                uzap1=uzap1+wpair(ii,iset)*
 325      &            harmonic(qpair(ii),qinpair(ii,iset))
 326             enddo
 327             ducartan(j,i)=(uzap2-uzap1)/(delta)
 328          enddo
 329       enddo
 330 c Calculating numerical gradients for dU/ddx
 331       do i=0,nres-1
 332          do j=1,3
 333            duxcartan(j,i)=0.0d0
 334          enddo
 335          do j=1,3
 336             cdummy(j,i)=dc(j,i+nres)
 337             dc(j,i+nres)=dc(j,i+nres)+delta
 338             call chainbuild_cart
 339             uzap2=0.0d0
 340             do ii=1,nfrag
 341                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 342      &         .true.,idummy,idummy)
 343                uzap2=uzap2+wfrag(ii,iset)*
 344      &            harmonic(qfrag(ii),qinfrag(ii,iset))
 345             enddo
 346             do ii=1,npair
 347                kstart=ifrag(1,ipair(1,ii,iset),iset)
 348                kend=ifrag(2,ipair(1,ii,iset),iset)
 349                lstart=ifrag(1,ipair(2,ii,iset),iset)
 350                lend=ifrag(2,ipair(2,ii,iset),iset)
 351                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 352                uzap2=uzap2+wpair(ii,iset)*
 353      &             harmonic(qpair(ii),qinpair(ii,iset))
 354             enddo
 355             dc(j,i+nres)=cdummy(j,i)
 356             call chainbuild_cart
 357             uzap1=0.0d0
 358              do ii=1,nfrag
 359                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 360      &         .true.,idummy,idummy)
 361                uzap1=uzap1+wfrag(ii,iset)*
 362      &            harmonic(qfrag(ii),qinfrag(ii,iset))
 363             enddo
 364             do ii=1,npair
 365                kstart=ifrag(1,ipair(1,ii,iset),iset)
 366                kend=ifrag(2,ipair(1,ii,iset),iset)
 367                lstart=ifrag(1,ipair(2,ii,iset),iset)
 368                lend=ifrag(2,ipair(2,ii,iset),iset)
 369                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 370                uzap1=uzap1+wpair(ii,iset)*
 371      &               harmonic(qpair(ii),qinpair(ii,iset))
 372             enddo
 373             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 374          enddo
 375       enddo
 376       write(iout,*) "Numerical dUconst/ddc backbone "
 377       do ii=0,nres
 378         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 379       enddo
 380       write(iout,*) "Numerical dUconst/ddx side-chain "
 381       do ii=1,nres
 382          write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 383       enddo
 384       return
 385       end
 386 c---------------------------------------------------------------------------
 387       double precision function qcontrib(il,jl,il1,jl1)
 388       implicit none
 389       include 'DIMENSIONS'
 390       include 'COMMON.IOUNITS'
 391       include 'COMMON.CHAIN'
 392       include 'COMMON.INTERACT'
 393       include 'COMMON.MD'
 394       integer i,j,k,il,jl,il1,jl1,nd
 395       double precision dist
 396       external dist
 397       double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
 398      &  fac1,ddave,ssij,ddqij
 399       logical lprn /.false./
 400       d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+
 401      &           (cref(2,jl)-cref(2,il))**2+
 402      &           (cref(3,jl)-cref(3,il))**2)
 403       dij1=dist(il,jl)
 404       ddave=(dij1-d0ij1)**2
 405       nd=1
 406       if (jl1.ne.jl) then
 407         d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+
 408      &           (cref(2,jl1)-cref(2,il))**2+
 409      &           (cref(3,jl1)-cref(3,il))**2)
 410         dij2=dist(il,jl1)
 411         ddave=ddave+(dij2-d0ij2)**2
 412         nd=nd+1
 413       endif
 414       if (il1.ne.il) then
 415         d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+
 416      &           (cref(2,jl)-cref(2,il1))**2+
 417      &           (cref(3,jl)-cref(3,il1))**2)
 418         dij3=dist(il1,jl)
 419         ddave=ddave+(dij3-d0ij3)**2
 420         nd=nd+1
 421       endif
 422       if (il1.ne.il .and. jl1.ne.jl) then
 423         d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+
 424      &           (cref(2,jl1)-cref(2,il1))**2+
 425      &           (cref(3,jl1)-cref(3,il1))**2)
 426         dij4=dist(il1,jl1)
 427         ddave=ddave+(dij4-d0ij4)**2
 428         nd=nd+1
 429       endif
 430       ddave=ddave/nd
 431       if (lprn) then
 432         write (iout,*) "il",il," jl",jl,
 433      &  " itype",itype(il),itype(jl)," nd",nd
 434         write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
 435      &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
 436         call flush(iout)
 437       endif
 438 c      ssij = (0.25d0*d0ij1)**2
 439       if (il.ne.il1 .and. jl.ne.jl1) then
 440         ssij = 16.0d0/(d0ij1*d0ij4)
 441       else
 442         ssij = 16.0d0/(d0ij1*d0ij1)
 443       endif
 444       qcontrib = dexp(-0.5d0*ddave*ssij)
 445 c Compute gradient
 446       fac1 = qcontrib*ssij/nd
 447       fac = fac1*(dij1-d0ij1)/dij1
 448       do k=1,3
 449         ddqij = (c(k,il)-c(k,jl))*fac
 450         dqwol(k,il)=dqwol(k,il)+ddqij
 451         dqwol(k,jl)=dqwol(k,jl)-ddqij
 452       enddo
 453       if (jl1.ne.jl) then
 454         fac = fac1*(dij2-d0ij2)/dij2
 455         do k=1,3
 456           ddqij = (c(k,il)-c(k,jl1))*fac
 457           dqwol(k,il)=dqwol(k,il)+ddqij
 458           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 459         enddo
 460       endif
 461       if (il1.ne.il) then
 462         fac = fac1*(dij3-d0ij3)/dij3
 463         do k=1,3
 464           ddqij = (c(k,il1)-c(k,jl))*fac
 465           dxqwol(k,il)=dxqwol(k,il)+ddqij
 466           dqwol(k,jl)=dqwol(k,jl)-ddqij
 467         enddo
 468       endif
 469       if (il1.ne.il .and. jl1.ne.jl) then
 470         fac = fac1*(dij4-d0ij4)/dij4
 471         do k=1,3
 472           ddqij = (c(k,il1)-c(k,jl1))*fac
 473           dxqwol(k,il)=dxqwol(k,il)+ddqij
 474           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 475         enddo
 476       endif
 477       return
 478       end