source/unres/src-HCD-5D/muca_md.F

   1       subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'COMMON.MUCA'
   5       include 'COMMON.MD'
   6 #ifdef LANG0
   7 #ifdef FIVEDIAG
   8       include 'COMMON.LANGEVIN.lang0.5diag'
   9 #else
  10       include 'COMMON.LANGEVIN.lang0'
  11 #endif
  12 #else
  13       include 'COMMON.LANGEVIN'
  14 #endif
  15       double precision remd_t_bath(maxprocs)
  16       double precision remd_ene(maxprocs)
  17       double precision muca_ene
  18       double precision betai,betaiex,delta
  19       integer i,iex
  20
  21       betai=1.0/(Rb*remd_t_bath(i))
  22       betaiex=1.0/(Rb*remd_t_bath(iex))
  23
  24       delta=betai*(muca_ene(remd_ene(iex),i,remd_t_bath)-
  25      &                muca_ene(remd_ene(i),i,remd_t_bath))
  26      &          -betaiex*(muca_ene(remd_ene(iex),iex,remd_t_bath)-
  27      &                muca_ene(remd_ene(i),iex,remd_t_bath))
  28
  29       return
  30       end
  31
  32       double precision function muca_ene(energy,i,remd_t_bath)
  33       implicit none
  34       include 'DIMENSIONS'
  35       include 'COMMON.MUCA'
  36       include 'COMMON.MD'
  37 #ifdef LANG0
  38 #ifdef FIVEDIAG
  39       include 'COMMON.LANGEVIN.lang0.5diag'
  40 #else
  41       include 'COMMON.LANGEVIN.lang0'
  42 #endif
  43 #else
  44       include 'COMMON.LANGEVIN'
  45 #endif
  46       double precision y,yp,energy
  47       double precision remd_t_bath(maxprocs)
  48       integer i
  49
  50       if (energy.lt.elowi(i)) then
  51         call splint(emuca,nemuca,nemuca2,nmuca,elowi(i),y,yp)
  52         muca_ene=remd_t_bath(i)*Rb*(yp*(energy-elowi(i))+y)
  53       elseif (energy.gt.ehighi(i)) then
  54         call splint(emuca,nemuca,nemuca2,nmuca,ehighi(i),y,yp)
  55         muca_ene=remd_t_bath(i)*Rb*(yp*(energy-ehighi(i))+y)
  56       else
  57         call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
  58         muca_ene=remd_t_bath(i)*Rb*y
  59       endif
  60       return
  61       end
  62
  63       subroutine read_muca
  64       implicit none
  65       include 'DIMENSIONS'
  66       include 'COMMON.MUCA'
  67       include 'COMMON.CONTROL'
  68       include 'COMMON.MD'
  69       include 'COMMON.REMD'
  70       include 'COMMON.SETUP'
  71       include 'COMMON.IOUNITS'
  72       double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
  73       integer i,j,k
  74       imtime=0
  75       do i=1,4*maxres
  76         hist(i)=0
  77       enddo
  78       if (modecalc.eq.14.and..not.remd_tlist) then
  79                 print *,"MUCAREMD works only with TLIST"
  80                 stop
  81       endif
  82       open(89,file='muca.input')
  83       read(89,*)
  84       read(89,*)
  85       if (modecalc.eq.14) then
  86         read(89,*) (elowi(i),ehighi(i),i=1,nrep)
  87        if (remd_mlist) then
  88         k=0
  89         do i=1,nrep
  90          do j=1,remd_m(i)
  91           i2rep(k)=i
  92           k=k+1
  93          enddo
  94         enddo
  95         elow=elowi(i2rep(me))
  96         ehigh=ehighi(i2rep(me))
  97         elowi(me+1)=elow
  98         ehighi(me+1)=ehigh
  99        else
 100         elow=elowi(me+1)
 101         ehigh=ehighi(me+1)
 102        endif
 103       else
 104         read(89,*) elow,ehigh
 105         elowi(1)=elow
 106         ehighi(1)=ehigh
 107       endif
 108       i=0
 109       do while(.true.)
 110        i=i+1
 111        read(89,*,end=100) emuca(i),nemuca(i)
 112 cd      nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
 113       enddo
 114  100  continue
 115       nmuca=i-1
 116       hbin=emuca(nmuca)-emuca(nmuca-1)
 117       write (iout,*) 'hbin',hbin
 118       write (iout,*) me,'elow,ehigh',elow,ehigh
 119       yp1=0
 120       ypn=0
 121       call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
 122       factor_min=0.0d0
 123       factor_min=muca_factor(ehigh)
 124       call print_muca
 125       return
 126       end
 127
 128
 129       subroutine print_muca
 130       implicit none
 131       include 'DIMENSIONS'
 132       include 'COMMON.MUCA'
 133       include 'COMMON.CONTROL'
 134       include 'COMMON.MD'
 135 #ifdef LANG0
 136 #ifdef FIVEDIAG
 137       include 'COMMON.LANGEVIN.lang0.5diag'
 138 #else
 139       include 'COMMON.LANGEVIN.lang0'
 140 #endif
 141 #else
 142       include 'COMMON.LANGEVIN'
 143 #endif
 144       include 'COMMON.REMD'
 145       include 'COMMON.SETUP'
 146       include 'COMMON.IOUNITS'
 147       double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
 148       double precision dummy(maxprocs)
 149       integer i,j,k
 150       if (remd_mlist) then
 151            k=0
 152            do i=1,nrep
 153             do j=1,remd_m(i)
 154              i2rep(k)=i
 155              k=k+1
 156             enddo
 157            enddo
 158       endif
 159
 160       do i=1,nmuca
 161 c       print *,'nemuca ',emuca(i),nemuca(i)
 162        do j=0,4
 163         x=emuca(i)+hbin/5*j
 164         if (modecalc.eq.14) then
 165          if (remd_mlist) then
 166           yp=muca_factor(x)*remd_t(i2rep(me))*Rb
 167           dummy(me+1)=remd_t(i2rep(me))
 168           y=muca_ene(x,me+1,dummy)
 169          else
 170           yp=muca_factor(x)*remd_t(me+1)*Rb
 171           y=muca_ene(x,me+1,remd_t)
 172          endif
 173          write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
 174      &      'muca factor ',x,yp,' muca ene',y
 175         else
 176          yp=muca_factor(x)*t_bath*Rb
 177          dummy(1)=t_bath
 178          y=muca_ene(x,1,dummy)
 179          write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
 180      &      'muca factor ',x,yp,' muca ene',y
 181         endif
 182        enddo
 183       enddo
 184       if(mucadyn.gt.0) then
 185        do i=1,nmuca
 186          write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
 187      &             imtime,emuca(i),nemuca(i)
 188        enddo
 189       endif
 190       return
 191       end
 192
 193       subroutine muca_update(energy)
 194       implicit none
 195       include 'DIMENSIONS'
 196       include 'COMMON.MUCA'
 197       include 'COMMON.CONTROL'
 198       include 'COMMON.MD'
 199       include 'COMMON.REMD'
 200       include 'COMMON.SETUP'
 201       include 'COMMON.IOUNITS'
 202       double precision energy
 203       double precision yp1,ypn
 204       integer i,j,k,ismooth,ist,ien
 205       logical lnotend
 206
 207       k=int((energy-emuca(1))/hbin)+1
 208
 209       IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
 210        if(energy.ge.ehigh)
 211      &        write (iout,*) 'MUCA reject',energy,emuca(k)
 212        if(energy.ge.ehigh.and.(energy-ehigh).lt.hbin) then
 213          write (iout,*) 'MUCA ehigh',energy,emuca(k)
 214          do i=k,nmuca
 215            hist(i)=hist(i)+1
 216          enddo
 217        endif
 218        if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1
 219       ELSE
 220        if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1
 221       ENDIF
 222       if(mod(imtime,mucadyn).eq.0) then
 223
 224          do i=1,nmuca
 225           IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
 226            nemuca(i)=nemuca(i)+dlog(hist(i)+1)
 227           ELSE
 228            if (hist(i).gt.0) hist(i)=dlog(hist(i))
 229            nemuca(i)=nemuca(i)+hist(i)
 230           ENDIF
 231           hist(i)=0
 232           write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
 233      &          imtime,emuca(i),nemuca(i)
 234          enddo
 235
 236
 237          lnotend=.true.
 238          ismooth=0
 239          ist=2
 240          ien=nmuca-1
 241         IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
 242 c         lnotend=.false.
 243 c         do i=1,nmuca-1
 244 c           do j=i+1,nmuca
 245 c            if(nemuca(j).lt.nemuca(i)) lnotend=.true.
 246 c           enddo
 247 c         enddo
 248          do while(lnotend)
 249           ismooth=ismooth+1
 250           write (iout,*) 'MUCA update smoothing',ist,ien
 251           do i=ist,ien
 252            nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
 253           enddo
 254           lnotend=.false.
 255           ist=0
 256           ien=0
 257           do i=1,nmuca-1
 258            do j=i+1,nmuca
 259             if(nemuca(j).lt.nemuca(i)) then
 260               lnotend=.true.
 261               if(ist.eq.0) ist=i-1
 262               if(ien.lt.j+1) ien=j+1
 263             endif
 264            enddo
 265           enddo
 266          enddo
 267         ENDIF
 268
 269          write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
 270          yp1=0
 271          ypn=0
 272          call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
 273          call print_muca
 274
 275       endif
 276       return
 277       end
 278
 279       double precision function muca_factor(energy)
 280       implicit none
 281       include 'DIMENSIONS'
 282       include 'COMMON.MUCA'
 283       double precision y,yp,energy
 284
 285       if (energy.lt.elow) then
 286         call splint(emuca,nemuca,nemuca2,nmuca,elow,y,yp)
 287       elseif (energy.gt.ehigh) then
 288         call splint(emuca,nemuca,nemuca2,nmuca,ehigh,y,yp)
 289       else
 290         call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
 291       endif
 292
 293       if(yp.ge.factor_min) then
 294        muca_factor=yp
 295       else
 296        muca_factor=factor_min
 297       endif
 298 cd      print *,'energy, muca_factor',energy,muca_factor
 299       return
 300       end
 301
 302
 303       SUBROUTINE spline(x,y,n,yp1,ypn,y2)
 304       implicit none
 305       INTEGER n,NMAX
 306       REAL*8 yp1,ypn,x(n),y(n),y2(n)
 307       PARAMETER (NMAX=500)
 308       INTEGER i,k
 309       REAL*8 p,qn,sig,un,u(NMAX)
 310       if (yp1.gt..99e30) then
 311       y2(1)=0.
 312       u(1)=0.
 313       else
 314          y2(1)=-0.5
 315          u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
 316       endif
 317       do i=2,n-1
 318          sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
 319          p=sig*y2(i-1)+2.
 320          y2(i)=(sig-1.)/p
 321          u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
 322      *        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
 323       enddo
 324       if (ypn.gt..99e30) then
 325          qn=0.
 326          un=0.
 327       else
 328          qn=0.5
 329          un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
 330       endif
 331       y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
 332       do k=n-1,1,-1
 333          y2(k)=y2(k)*y2(k+1)+u(k)
 334       enddo
 335       return
 336       END
 337
 338
 339       SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
 340       implicit none
 341       INTEGER n
 342       REAL*8 x,y,xa(n),y2a(n),ya(n),yp
 343       INTEGER k,khi,klo
 344       REAL*8 a,b,h
 345       klo=1
 346       khi=n
 347  1    if (khi-klo.gt.1) then
 348          k=(khi+klo)/2
 349          if (xa(k).gt.x) then
 350             khi=k
 351          else
 352             klo=k
 353          endif
 354          goto 1
 355       endif
 356       h=xa(khi)-xa(klo)
 357       if (h.eq.0.) pause 'bad xa input in splint'
 358       a=(xa(khi)-x)/h
 359       b=(x-xa(klo))/h
 360       y=a*ya(klo)+b*ya(khi)+
 361      *     ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
 362       yp=-ya(klo)/h+ya(khi)/h-3*(a**2)*y2a(klo)*h/6.
 363      +     +(3*(b**2)-1)*y2a(khi)*h/6.
 364       return
 365       END