source/unres/src_MD/src/q_measure3.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,iset),ifrag(2,i,iset),.true.,idummy,idummy)
 199       enddo
 200 c      stop
 201       do i=1,npair
 202          kstart=ifrag(1,ipair(1,i,iset),iset)
 203          kend=ifrag(2,ipair(1,i,iset),iset)
 204          lstart=ifrag(1,ipair(2,i,iset),iset)
 205          lend=ifrag(2,ipair(2,i,iset),iset)
 206          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
 207          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
 208 c  Calculating dU/dQ
 209          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
 210 c Calculating dQ/dXi
 211          do ii=0,nres
 212             do j=1,3
 213                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 214                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 215             enddo
 216          enddo
 217       enddo
 218 c      write(iout,*) "Uconst inside subroutine ", Uconst
 219 c Transforming the gradients from Cs to dCs for the backbone
 220       do i=0,nres
 221          do j=i+1,nres
 222            do k=1,3
 223              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 224            enddo
 225          enddo
 226       enddo
 227 c  Transforming the gradients from Cs to dCs for the side chains
 228       do i=1,nres
 229          do j=1,3
 230            dudxconst(j,i)=duxconst(j,i)
 231          enddo
 232       enddo
 233 c      write(iout,*) "dU/dc backbone "
 234 c       do ii=0,nres
 235 c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
 236 c      enddo
 237 c      write(iout,*) "dU/dX side chain "
 238 c      do ii=1,nres
 239 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 240 c      enddo
 241 c      write(iout,*) "dU/ddc backbone "
 242 c       do ii=0,nres
 243 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 244 c      enddo
 245 c      write(iout,*) "dU/ddX side chain "
 246 c      do ii=1,nres
 247 c            write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
 248 c      enddo
 249 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 250 c      call dEconstrQ_num
 251       return
 252       end
 253 c-----------------------------------------------------------------------
 254       subroutine dEconstrQ_num
 255 c Calculating numerical dUconst/ddc and dUconst/ddx
 256       implicit real*8 (a-h,o-z)
 257       include 'DIMENSIONS'
 258       include 'COMMON.CONTROL'
 259       include 'COMMON.VAR'
 260       include 'COMMON.MD'
 261 #ifndef LANG0
 262       include 'COMMON.LANGEVIN'
 263 #else
 264       include 'COMMON.LANGEVIN.lang0'
 265 #endif
 266       include 'COMMON.CHAIN'
 267       include 'COMMON.DERIV'
 268       include 'COMMON.GEO'
 269       include 'COMMON.LOCAL'
 270       include 'COMMON.INTERACT'
 271       include 'COMMON.IOUNITS'
 272       include 'COMMON.NAMES'
 273       include 'COMMON.TIME1'
 274       double precision uzap1,uzap2
 275       double precision dUcartan(3,0:MAXRES)
 276      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
 277       integer kstart,kend,lstart,lend,idummy
 278       double precision delta /1.0d-7/
 279 c     For the backbone
 280       do i=0,nres-1
 281          do j=1,3
 282             dUcartan(j,i)=0.0d0
 283             cdummy(j,i)=dc(j,i)
 284             dc(j,i)=dc(j,i)+delta
 285             call chainbuild_cart
 286             uzap2=0.0d0
 287             do ii=1,nfrag
 288                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 289      &           .true.,idummy,idummy)
 290                uzap2=uzap2+wfrag(ii,iset)*
 291      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 292             enddo
 293             do ii=1,npair
 294                kstart=ifrag(1,ipair(1,ii,iset),iset)
 295                kend=ifrag(2,ipair(1,ii,iset),iset)
 296                lstart=ifrag(1,ipair(2,ii,iset),iset)
 297                lend=ifrag(2,ipair(2,ii,iset),iset)
 298                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 299                uzap2=uzap2+wpair(ii,iset)*
 300      &                harmonic(qpair(ii),qinpair(ii,iset))
 301             enddo
 302             dc(j,i)=cdummy(j,i)
 303             call chainbuild_cart
 304             uzap1=0.0d0
 305              do ii=1,nfrag
 306                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 307      &           .true.,idummy,idummy)
 308                uzap1=uzap1+wfrag(ii,iset)*
 309      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 310             enddo
 311             do ii=1,npair
 312                kstart=ifrag(1,ipair(1,ii,iset),iset)
 313                kend=ifrag(2,ipair(1,ii,iset),iset)
 314                lstart=ifrag(1,ipair(2,ii,iset),iset)
 315                lend=ifrag(2,ipair(2,ii,iset),iset)
 316                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 317                uzap1=uzap1+wpair(ii,iset)*
 318      &                harmonic(qpair(ii),qinpair(ii,iset))
 319             enddo
 320             ducartan(j,i)=(uzap2-uzap1)/(delta)
 321          enddo
 322       enddo
 323 c Calculating numerical gradients for dU/ddx
 324       do i=0,nres-1
 325          do j=1,3
 326            duxcartan(j,i)=0.0d0
 327          enddo
 328          do j=1,3
 329             cdummy(j,i)=dc(j,i+nres)
 330             dc(j,i+nres)=dc(j,i+nres)+delta
 331             call chainbuild_cart
 332             uzap2=0.0d0
 333             do ii=1,nfrag
 334                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 335      &           .true.,idummy,idummy)
 336                uzap2=uzap2+wfrag(ii,iset)*
 337      &                harmonic(qfrag(ii),qinfrag(ii,iset))
 338             enddo
 339             do ii=1,npair
 340                kstart=ifrag(1,ipair(1,ii,iset),iset)
 341                kend=ifrag(2,ipair(1,ii,iset),iset)
 342                lstart=ifrag(1,ipair(2,ii,iset),iset)
 343                lend=ifrag(2,ipair(2,ii,iset),iset)
 344                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 345                uzap2=uzap2+wpair(ii,iset)*
 346      &                harmonic(qpair(ii),qinpair(ii,iset))
 347             enddo
 348             dc(j,i+nres)=cdummy(j,i)
 349             call chainbuild_cart
 350             uzap1=0.0d0
 351              do ii=1,nfrag
 352                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
 353      &           .true.,idummy,idummy)
 354                uzap1=uzap1+wfrag(ii,iset)*
 355      &                 harmonic(qfrag(ii),qinfrag(ii,iset))
 356             enddo
 357             do ii=1,npair
 358                kstart=ifrag(1,ipair(1,ii,iset),iset)
 359                kend=ifrag(2,ipair(1,ii,iset),iset)
 360                lstart=ifrag(1,ipair(2,ii,iset),iset)
 361                lend=ifrag(2,ipair(2,ii,iset),iset)
 362                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 363                uzap1=uzap1+wpair(ii,iset)*
 364      &                harmonic(qpair(ii),qinpair(ii,iset))
 365             enddo
 366             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 367          enddo
 368       enddo
 369       write(iout,*) "Numerical dUconst/ddc backbone "
 370       do ii=0,nres
 371         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 372       enddo
 373       write(iout,*) "Numerical dUconst/ddx side-chain "
 374       do ii=1,nres
 375          write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 376       enddo
 377       return
 378       end
 379 c---------------------------------------------------------------------------
 380       double precision function qcontrib(il,jl,il1,jl1)
 381       implicit none
 382       include 'DIMENSIONS'
 383       include 'COMMON.IOUNITS'
 384       include 'COMMON.CHAIN'
 385       include 'COMMON.INTERACT'
 386       include 'COMMON.MD'
 387       include 'COMMON.LOCAL'
 388       integer i,j,k,il,jl,il1,jl1,nd,itl,jtl
 389       double precision dist
 390       external dist
 391       double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120
 392      &  ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12
 393       double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3),
 394      &  aux1,aux2
 395       double precision scalar
 396       external scalar
 397       logical lprn /.false./
 398       if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1
 399       d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 400      &           (cref(2,jl)-cref(2,il))**2+
 401      &           (cref(3,jl)-cref(3,il))**2)
 402       dij=dist(il,jl)
 403       dij1=dist(il1,jl1)
 404       do i=1,3
 405         er(i)=(c(i,jl1)-c(i,il1))/dij1
 406       enddo
 407       do i=1,3
 408         er0(i)=cref(i,jl1)-cref(i,il1)
 409       enddo
 410       d0ij1=dsqrt(scalar(er0,er0))
 411       do i=1,3
 412         er0(i)=er0(i)/d0ij1
 413       enddo
 414       if (il.ne.il1 .or. jl.ne.jl1) then
 415         ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2)
 416         nd=2
 417       else
 418         ddave=(dij-d0ij)**2
 419         nd=1
 420       endif
 421       if (il.ne.il1) then
 422         do i=1,3
 423           u(i)=cref(i,il1)-cref(i,il)
 424         enddo
 425         d0ii1=dsqrt(scalar(u,u))
 426         do i=1,3
 427           u(i)=u(i)/d0ii1
 428         enddo
 429         if (lprn) then
 430         write (iout,*) "u",(u(i),i=1,3)
 431         write (iout,*) "er0",(er0(i),i=1,3)
 432         om10=scalar(er0,u)
 433         om1=scalar(er,dc_norm(1,il1))
 434         write (iout,*) "om10",om10," om1",om1
 435         endif
 436       else
 437         om1=0.0d0
 438         om10=0.0d0
 439       endif
 440       if (jl.ne.jl1) then
 441         do i=1,3
 442           v(i)=cref(i,jl1)-cref(i,jl)
 443         enddo
 444         d0jj1=dsqrt(scalar(v,v))
 445         do i=1,3
 446           v(i)=v(i)/d0jj1
 447         enddo
 448         if (lprn) then
 449         write (iout,*) "v",(v(i),i=1,3)
 450         write (iout,*) "er0",(er0(i),i=1,3)
 451         om20=scalar(er,v)
 452         om2=scalar(er,dc_norm(1,jl1))
 453         write (iout,*) "om20",om20," om2",om2
 454         endif
 455       else
 456         om2=0.0d0
 457         om20=0.0d0
 458       endif
 459       if (il.ne.il1 .and. jl.ne.jl1) then
 460         om120=scalar(u,v)
 461         om12=scalar(dc_norm(1,il1),dc_norm(1,jl1))
 462       else
 463         om12=0.0d0
 464         om120=0.0d0
 465       endif
 466       if (lprn) then
 467         write (iout,*) "il",il," jl",jl,itype(il),itype(jl)
 468         write (iout,*)"d0ij",d0ij," om10",om10," om20",om20,
 469      &   " om120",om120,
 470      &  " dij",dij," om1",om1," om2",om2," om12",om12
 471         call flush(iout)
 472       endif
 473       ssij = 16.0d0/(d0ij*d0ij)
 474       qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2
 475      &                       +(om2-om20)**2+(om12-om120)**2)))
 476       if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib
 477 c      qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2)
 478 c      qcontrib = dexp(-0.5d0*(ddave*ssij))
 479 c Compute gradient - radial component
 480       fac1 = qcontrib*ssij/nd
 481       fac = fac1*(dij-d0ij)/dij
 482       do k=1,3
 483         ddqij = (c(k,il)-c(k,jl))*fac
 484         dqwol(k,il)=dqwol(k,il)+ddqij
 485         dqwol(k,jl)=dqwol(k,jl)-ddqij
 486       enddo
 487       if (il1.ne.il .or. jl1.ne.jl) then
 488         fac = fac1*(dij1-d0ij1)/dij1
 489         do k=1,3
 490           ddqij = (c(k,il1)-c(k,jl1))*fac
 491           if (il1.ne.il) then
 492             dxqwol(k,il)=dxqwol(k,il)+ddqij
 493           else
 494             dqwol(k,il)=dqwol(k,il)+ddqij
 495           endif
 496           if (jl1.ne.jl) then
 497             dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 498           else
 499             dqwol(k,jl)=dqwol(k,jl)-ddqij
 500           endif
 501         enddo
 502       endif
 503 c      return
 504 c Orientational contributions
 505       rij=1.0d0/dij1
 506       eom1=qcontrib*(om1-om10)
 507       eom2=qcontrib*(om2-om20)
 508       eom12=qcontrib*(om12-om120)
 509       do k=1,3
 510         dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k))
 511         dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k))
 512       enddo
 513       do k=1,3
 514         ddqij=eom1*dcosom1(k)+eom2*dcosom2(k)
 515         aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
 516      &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
 517         aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
 518      &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
 519         dqwol(k,il)=dqwol(k,il)-ddqij-aux1
 520         dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2
 521         dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1
 522 c     &            +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
 523 c     &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
 524         dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2
 525 c     &            +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
 526 c     &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
 527       enddo
 528       return
 529       end