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