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