source/unres/src_MD-M/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       include 'COMMON.LANGEVIN.lang0'
 160 #endif
 161       include 'COMMON.CHAIN'
 162       include 'COMMON.DERIV'
 163       include 'COMMON.GEO'
 164       include 'COMMON.LOCAL'
 165       include 'COMMON.INTERACT'
 166       include 'COMMON.IOUNITS'
 167       include 'COMMON.NAMES'
 168       include 'COMMON.TIME1'
 169       double precision uzap1,uzap2,hm1,hm2,hmnum
 170       double precision ucdelan,dUcartan(3,0:MAXRES)
 171      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
 172      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
 173       integer kstart,kend,lstart,lend,idummy
 174       double precision delta /1.0d-7/
 175       do i=0,nres
 176          do j=1,3
 177             duconst(j,i)=0.0d0
 178             dudconst(j,i)=0.0d0
 179             duxconst(j,i)=0.0d0
 180             dudxconst(j,i)=0.0d0
 181          enddo
 182       enddo
 183       Uconst=0.0d0
 184       do i=1,nfrag
 185          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 186      &   ,idummy,idummy)
 187          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
 188 c Calculating the derivatives of Constraint energy with respect to Q
 189          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
 190 c Calculating the derivatives of Q with respect to cartesian coordinates
 191          do ii=0,nres
 192             do j=1,3
 193                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
 194                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
 195             enddo
 196          enddo
 197 c      write (iout,*) "Calling qwol_num"
 198 c      call qwol_num(ifrag(1,i),ifrag(2,i),.true.,idummy,idummy)
 199       enddo
 200       do i=1,npair
 201          kstart=ifrag(1,ipair(1,i,iset),iset)
 202          kend=ifrag(2,ipair(1,i,iset),iset)
 203          lstart=ifrag(1,ipair(2,i,iset),iset)
 204          lend=ifrag(2,ipair(2,i,iset),iset)
 205          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
 206          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
 207 c  Calculating dU/dQ
 208          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
 209 c Calculating dQ/dXi
 210          do ii=0,nres
 211             do j=1,3
 212                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 213                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 214             enddo
 215          enddo
 216       enddo
 217 c      write(iout,*) "Uconst inside subroutine ", Uconst
 218 c Transforming the gradients from Cs to dCs for the backbone
 219       do i=0,nres
 220          do j=i+1,nres
 221            do k=1,3
 222              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 223            enddo
 224          enddo
 225       enddo
 226 c  Transforming the gradients from Cs to dCs for the side chains
 227       do i=1,nres
 228          do j=1,3
 229            dudxconst(j,i)=duxconst(j,i)
 230          enddo
 231       enddo
 232 c      write(iout,*) "dU/dc backbone "
 233 c       do ii=0,nres
 234 c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
 235 c      enddo
 236 c      write(iout,*) "dU/dX side chain "
 237 c      do ii=1,nres
 238 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 239 c      enddo
 240 c      write(iout,*) "dU/ddc backbone "
 241 c       do ii=0,nres
 242 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 243 c      enddo
 244 c      write(iout,*) "dU/ddX side chain "
 245 c      do ii=1,nres
 246 c            write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
 247 c      enddo
 248 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 249 c      call dEconstrQ_num
 250       return
 251       end
 252 c-----------------------------------------------------------------------
 253       subroutine dEconstrQ_num
 254 c Calculating numerical dUconst/ddc and dUconst/ddx
 255       implicit real*8 (a-h,o-z)
 256       include 'DIMENSIONS'
 257       include 'COMMON.CONTROL'
 258       include 'COMMON.VAR'
 259       include 'COMMON.MD'
 260 #ifndef LANG0
 261       include 'COMMON.LANGEVIN'
 262 #else
 263       include 'COMMON.LANGEVIN.lang0'
 264 #endif
 265       include 'COMMON.CHAIN'
 266       include 'COMMON.DERIV'
 267       include 'COMMON.GEO'
 268       include 'COMMON.LOCAL'
 269       include 'COMMON.INTERACT'
 270       include 'COMMON.IOUNITS'
 271       include 'COMMON.NAMES'
 272       include 'COMMON.TIME1'
 273       double precision uzap1,uzap2
 274       double precision dUcartan(3,0:MAXRES)
 275      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
 276       integer kstart,kend,lstart,lend,idummy
 277       double precision delta /1.0d-7/
 278 c     For the backbone
 279       do i=0,nres-1
 280          do j=1,3
 281             dUcartan(j,i)=0.0d0
 282             cdummy(j,i)=dc(j,i)
 283             dc(j,i)=dc(j,i)+delta
 284             call chainbuild_cart
 285             uzap2=0.0d0
 286             do ii=1,nfrag
 287                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 288      &         .true.,idummy,idummy)
 289                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
 290      &            qinfrag(ii,iset))
 291             enddo
 292             do ii=1,npair
 293                kstart=ifrag(1,ipair(1,ii,iset),iset)
 294                kend=ifrag(2,ipair(1,ii,iset),iset)
 295                lstart=ifrag(1,ipair(2,ii,iset),iset)
 296                lend=ifrag(2,ipair(2,ii,iset),iset)
 297                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 298                uzap2=uzap2+wpair(ii,iset)*
 299      &             harmonic(qpair(ii),qinpair(ii,iset))
 300             enddo
 301             dc(j,i)=cdummy(j,i)
 302             call chainbuild_cart
 303             uzap1=0.0d0
 304              do ii=1,nfrag
 305                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 306      &         .true.,idummy,idummy)
 307                uzap1=uzap1+wfrag(ii,iset)*
 308      &                 harmonic(qfrag(ii),qinfrag(ii,iset))
 309             enddo
 310             do ii=1,npair
 311                kstart=ifrag(1,ipair(1,ii,iset),iset)
 312                kend=ifrag(2,ipair(1,ii,iset),iset)
 313                lstart=ifrag(1,ipair(2,ii,iset),iset)
 314                lend=ifrag(2,ipair(2,ii,iset),iset)
 315                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 316                uzap1=uzap1+wpair(ii,iset)*
 317      &            harmonic(qpair(ii),qinpair(ii,iset))
 318             enddo
 319             ducartan(j,i)=(uzap2-uzap1)/(delta)
 320          enddo
 321       enddo
 322 c Calculating numerical gradients for dU/ddx
 323       do i=0,nres-1
 324          do j=1,3
 325            duxcartan(j,i)=0.0d0
 326          enddo
 327          do j=1,3
 328             cdummy(j,i)=dc(j,i+nres)
 329             dc(j,i+nres)=dc(j,i+nres)+delta
 330             call chainbuild_cart
 331             uzap2=0.0d0
 332             do ii=1,nfrag
 333                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 334      &         .true.,idummy,idummy)
 335                uzap2=uzap2+wfrag(ii,iset)*
 336      &            harmonic(qfrag(ii),qinfrag(ii,iset))
 337             enddo
 338             do ii=1,npair
 339                kstart=ifrag(1,ipair(1,ii,iset),iset)
 340                kend=ifrag(2,ipair(1,ii,iset),iset)
 341                lstart=ifrag(1,ipair(2,ii,iset),iset)
 342                lend=ifrag(2,ipair(2,ii,iset),iset)
 343                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 344                uzap2=uzap2+wpair(ii,iset)*
 345      &             harmonic(qpair(ii),qinpair(ii,iset))
 346             enddo
 347             dc(j,i+nres)=cdummy(j,i)
 348             call chainbuild_cart
 349             uzap1=0.0d0
 350              do ii=1,nfrag
 351                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 352      &         .true.,idummy,idummy)
 353                uzap1=uzap1+wfrag(ii,iset)*
 354      &            harmonic(qfrag(ii),qinfrag(ii,iset))
 355             enddo
 356             do ii=1,npair
 357                kstart=ifrag(1,ipair(1,ii,iset),iset)
 358                kend=ifrag(2,ipair(1,ii,iset),iset)
 359                lstart=ifrag(1,ipair(2,ii,iset),iset)
 360                lend=ifrag(2,ipair(2,ii,iset),iset)
 361                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 362                uzap1=uzap1+wpair(ii,iset)*
 363      &               harmonic(qpair(ii),qinpair(ii,iset))
 364             enddo
 365             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 366          enddo
 367       enddo
 368       write(iout,*) "Numerical dUconst/ddc backbone "
 369       do ii=0,nres
 370         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 371       enddo
 372       write(iout,*) "Numerical dUconst/ddx side-chain "
 373       do ii=1,nres
 374          write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 375       enddo
 376       return
 377       end
 378 c---------------------------------------------------------------------------
 379       double precision function qcontrib(il,jl,il1,jl1)
 380       implicit none
 381       include 'DIMENSIONS'
 382       include 'COMMON.IOUNITS'
 383       include 'COMMON.CHAIN'
 384       include 'COMMON.INTERACT'
 385       include 'COMMON.MD'
 386       integer i,j,k,il,jl,il1,jl1,nd
 387       double precision dist
 388       external dist
 389       double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
 390      &  fac1,ddave,ssij,ddqij
 391       logical lprn /.false./
 392       d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+
 393      &           (cref(2,jl)-cref(2,il))**2+
 394      &           (cref(3,jl)-cref(3,il))**2)
 395       dij1=dist(il,jl)
 396       ddave=(dij1-d0ij1)**2
 397       nd=1
 398       if (jl1.ne.jl) then
 399         d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+
 400      &           (cref(2,jl1)-cref(2,il))**2+
 401      &           (cref(3,jl1)-cref(3,il))**2)
 402         dij2=dist(il,jl1)
 403         ddave=ddave+(dij2-d0ij2)**2
 404         nd=nd+1
 405       endif
 406       if (il1.ne.il) then
 407         d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+
 408      &           (cref(2,jl)-cref(2,il1))**2+
 409      &           (cref(3,jl)-cref(3,il1))**2)
 410         dij3=dist(il1,jl)
 411         ddave=ddave+(dij3-d0ij3)**2
 412         nd=nd+1
 413       endif
 414       if (il1.ne.il .and. jl1.ne.jl) then
 415         d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+
 416      &           (cref(2,jl1)-cref(2,il1))**2+
 417      &           (cref(3,jl1)-cref(3,il1))**2)
 418         dij4=dist(il1,jl1)
 419         ddave=ddave+(dij4-d0ij4)**2
 420         nd=nd+1
 421       endif
 422       ddave=ddave/nd
 423       if (lprn) then
 424         write (iout,*) "il",il," jl",jl,
 425      &  " itype",itype(il),itype(jl)," nd",nd
 426         write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
 427      &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
 428         call flush(iout)
 429       endif
 430 c      ssij = (0.25d0*d0ij1)**2
 431       if (il.ne.il1 .and. jl.ne.jl1) then
 432         ssij = 16.0d0/(d0ij1*d0ij4)
 433       else
 434         ssij = 16.0d0/(d0ij1*d0ij1)
 435       endif
 436       qcontrib = dexp(-0.5d0*ddave*ssij)
 437 c Compute gradient
 438       fac1 = qcontrib*ssij/nd
 439       fac = fac1*(dij1-d0ij1)/dij1
 440       do k=1,3
 441         ddqij = (c(k,il)-c(k,jl))*fac
 442         dqwol(k,il)=dqwol(k,il)+ddqij
 443         dqwol(k,jl)=dqwol(k,jl)-ddqij
 444       enddo
 445       if (jl1.ne.jl) then
 446         fac = fac1*(dij2-d0ij2)/dij2
 447         do k=1,3
 448           ddqij = (c(k,il)-c(k,jl1))*fac
 449           dqwol(k,il)=dqwol(k,il)+ddqij
 450           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 451         enddo
 452       endif
 453       if (il1.ne.il) then
 454         fac = fac1*(dij3-d0ij3)/dij3
 455         do k=1,3
 456           ddqij = (c(k,il1)-c(k,jl))*fac
 457           dxqwol(k,il)=dxqwol(k,il)+ddqij
 458           dqwol(k,jl)=dqwol(k,jl)-ddqij
 459         enddo
 460       endif
 461       if (il1.ne.il .and. jl1.ne.jl) then
 462         fac = fac1*(dij4-d0ij4)/dij4
 463         do k=1,3
 464           ddqij = (c(k,il1)-c(k,jl1))*fac
 465           dxqwol(k,il)=dxqwol(k,il)+ddqij
 466           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 467         enddo
 468       endif
 469       return
 470       end