source/unres/src-HCD-5D/q_measure-02.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.02d0*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.02d0*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 #ifdef FIVEDIAG
 244       include 'COMMON.LANGEVIN.lang0.5diag'
 245 #else
 246       include 'COMMON.LANGEVIN.lang0'
 247 #endif
 248 #endif
 249       include 'COMMON.CHAIN'
 250       include 'COMMON.DERIV'
 251       include 'COMMON.GEO'
 252       include 'COMMON.LOCAL'
 253       include 'COMMON.INTERACT'
 254       include 'COMMON.IOUNITS'
 255       include 'COMMON.NAMES'
 256       include 'COMMON.TIME1'
 257       double precision uzap1,uzap2,hm1,hm2,hmnum
 258       double precision ucdelan,dUcartan(3,0:MAXRES)
 259      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
 260      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
 261       integer kstart,kend,lstart,lend,idummy
 262       double precision delta /1.0d-7/
 263       do i=0,nres
 264          do j=1,3
 265             duconst(j,i)=0.0d0
 266             dudconst(j,i)=0.0d0
 267             duxconst(j,i)=0.0d0
 268             dudxconst(j,i)=0.0d0
 269          enddo
 270       enddo
 271       Uconst=0.0d0
 272       do i=1,nfrag
 273          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 274      &    ,idummy,idummy)
 275          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
 276 c Calculating the derivatives of Constraint energy with respect to Q
 277          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
 278      &     qinfrag(i,iset))
 279 c         hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
 280 c        hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
 281 c         hmnum=(hm2-hm1)/delta
 282 c         write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
 283 c     &   qinfrag(i,iset))
 284 c         write(iout,*) "harmonicnum frag", hmnum
 285 c Calculating the derivatives of Q with respect to cartesian coordinates
 286          call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 287      &   ,idummy,idummy)
 288 c         write(iout,*) "dqwol "
 289 c         do ii=1,nres
 290 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
 291 c         enddo
 292 c         write(iout,*) "dxqwol "
 293 c         do ii=1,nres
 294 c           write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
 295 c         enddo
 296 c Calculating numerical gradients of dU/dQi and dQi/dxi
 297 c        call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
 298 c     &  ,idummy,idummy)
 299 c  The gradients of Uconst in Cs
 300          do ii=0,nres
 301             do j=1,3
 302                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
 303                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
 304             enddo
 305          enddo
 306       enddo
 307       do i=1,npair
 308          kstart=ifrag(1,ipair(1,i,iset),iset)
 309          kend=ifrag(2,ipair(1,i,iset),iset)
 310          lstart=ifrag(1,ipair(2,i,iset),iset)
 311          lend=ifrag(2,ipair(2,i,iset),iset)
 312          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
 313          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
 314 c  Calculating dU/dQ
 315          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
 316 c         hm1=harmonic(qpair(i),qinpair(i,iset))
 317 c        hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
 318 c         hmnum=(hm2-hm1)/delta
 319 c         write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
 320 c     &   qinpair(i,iset))
 321 c         write(iout,*) "harmonicnum pair ", hmnum
 322 c Calculating dQ/dXi
 323          call qwolynes_prim(kstart,kend,.false.
 324      &   ,lstart,lend)
 325 c         write(iout,*) "dqwol "
 326 c         do ii=1,nres
 327 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
 328 c         enddo
 329 c         write(iout,*) "dxqwol "
 330 c         do ii=1,nres
 331 c          write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
 332 c        enddo
 333 c Calculating numerical gradients
 334 c        call qwol_num(kstart,kend,.false.
 335 c     &  ,lstart,lend)
 336 c The gradients of Uconst in Cs
 337          do ii=0,nres
 338             do j=1,3
 339                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
 340                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
 341             enddo
 342          enddo
 343       enddo
 344 c      write(iout,*) "Uconst inside subroutine ", Uconst
 345 c Transforming the gradients from Cs to dCs for the backbone
 346       do i=0,nres
 347          do j=i+1,nres
 348            do k=1,3
 349              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
 350            enddo
 351          enddo
 352       enddo
 353 c  Transforming the gradients from Cs to dCs for the side chains
 354       do i=1,nres
 355          do j=1,3
 356            dudxconst(j,i)=duxconst(j,i)
 357          enddo
 358       enddo
 359 c      write(iout,*) "dU/ddc backbone "
 360 c       do ii=0,nres
 361 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
 362 c      enddo
 363 c      write(iout,*) "dU/ddX side chain "
 364 c      do ii=1,nres
 365 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
 366 c      enddo
 367 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
 368 c      call dEconstrQ_num
 369       return
 370       end
 371 c-----------------------------------------------------------------------
 372       subroutine dEconstrQ_num
 373 c Calculating numerical dUconst/ddc and dUconst/ddx
 374       implicit real*8 (a-h,o-z)
 375       include 'DIMENSIONS'
 376       include 'COMMON.CONTROL'
 377       include 'COMMON.VAR'
 378       include 'COMMON.MD'
 379 #ifndef LANG0
 380       include 'COMMON.LANGEVIN'
 381 #else
 382 #ifdef FIVEDIAG
 383       include 'COMMON.LANGEVIN.lang0.5diag'
 384 #else
 385       include 'COMMON.LANGEVIN.lang0'
 386 #endif
 387 #endif
 388       include 'COMMON.CHAIN'
 389       include 'COMMON.DERIV'
 390       include 'COMMON.GEO'
 391       include 'COMMON.LOCAL'
 392       include 'COMMON.INTERACT'
 393       include 'COMMON.IOUNITS'
 394       include 'COMMON.NAMES'
 395       include 'COMMON.TIME1'
 396       double precision uzap1,uzap2
 397       double precision dUcartan(3,0:MAXRES)
 398      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
 399       integer kstart,kend,lstart,lend,idummy
 400       double precision delta /1.0d-7/
 401 c     For the backbone
 402       do i=0,nres-1
 403          do j=1,3
 404             dUcartan(j,i)=0.0d0
 405             cdummy(j,i)=dc(j,i)
 406             dc(j,i)=dc(j,i)+delta
 407             call chainbuild_cart
 408             uzap2=0.0d0
 409             do ii=1,nfrag
 410              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
 411      &         ,idummy,idummy)
 412                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
 413      &          qinfrag(ii,iset))
 414             enddo
 415             do ii=1,npair
 416                kstart=ifrag(1,ipair(1,ii,iset),iset)
 417                kend=ifrag(2,ipair(1,ii,iset),iset)
 418                lstart=ifrag(1,ipair(2,ii,iset),iset)
 419                lend=ifrag(2,ipair(2,ii,iset),iset)
 420                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 421                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
 422      &           qinpair(ii,iset))
 423             enddo
 424             dc(j,i)=cdummy(j,i)
 425             call chainbuild_cart
 426             uzap1=0.0d0
 427              do ii=1,nfrag
 428              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
 429      &         ,idummy,idummy)
 430                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
 431      &          qinfrag(ii,iset))
 432             enddo
 433             do ii=1,npair
 434                kstart=ifrag(1,ipair(1,ii,iset),iset)
 435                kend=ifrag(2,ipair(1,ii,iset),iset)
 436                lstart=ifrag(1,ipair(2,ii,iset),iset)
 437                lend=ifrag(2,ipair(2,ii,iset),iset)
 438                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 439                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
 440      &          qinpair(ii,iset))
 441             enddo
 442             ducartan(j,i)=(uzap2-uzap1)/(delta)
 443          enddo
 444       enddo
 445 c Calculating numerical gradients for dU/ddx
 446       do i=0,nres-1
 447          duxcartan(j,i)=0.0d0
 448          do j=1,3
 449             cdummy(j,i)=dc(j,i+nres)
 450             dc(j,i+nres)=dc(j,i+nres)+delta
 451             call chainbuild_cart
 452             uzap2=0.0d0
 453             do ii=1,nfrag
 454              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
 455      &         ,idummy,idummy)
 456                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
 457      &          qinfrag(ii,iset))
 458             enddo
 459             do ii=1,npair
 460                kstart=ifrag(1,ipair(1,ii,iset),iset)
 461                kend=ifrag(2,ipair(1,ii,iset),iset)
 462                lstart=ifrag(1,ipair(2,ii,iset),iset)
 463                lend=ifrag(2,ipair(2,ii,iset),iset)
 464                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 465                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
 466      &          qinpair(ii,iset))
 467             enddo
 468             dc(j,i+nres)=cdummy(j,i)
 469             call chainbuild_cart
 470             uzap1=0.0d0
 471              do ii=1,nfrag
 472                qfrag(ii)=qwolynes(ifrag(1,ii,iset),
 473      &          ifrag(2,ii,iset),.true.,idummy,idummy)
 474                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
 475      &          qinfrag(ii,iset))
 476             enddo
 477             do ii=1,npair
 478                kstart=ifrag(1,ipair(1,ii,iset),iset)
 479                kend=ifrag(2,ipair(1,ii,iset),iset)
 480                lstart=ifrag(1,ipair(2,ii,iset),iset)
 481                lend=ifrag(2,ipair(2,ii,iset),iset)
 482                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
 483                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
 484      &          qinpair(ii,iset))
 485             enddo
 486             duxcartan(j,i)=(uzap2-uzap1)/(delta)
 487          enddo
 488       enddo
 489       write(iout,*) "Numerical dUconst/ddc backbone "
 490       do ii=0,nres
 491         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
 492       enddo
 493 c      write(iout,*) "Numerical dUconst/ddx side-chain "
 494 c      do ii=1,nres
 495 c         write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
 496 c      enddo
 497       return
 498       end
 499 c---------------------------------------------------------------------------