source/unres/src_MD-M/muca_md.f

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