source/unres/muca_md.F90

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