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