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