--- /dev/null
+ subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ include 'COMMON.MD'
+ double precision remd_t_bath(maxprocs)
+ double precision remd_ene(maxprocs)
+ double precision muca_ene
+ double precision betai,betaiex,delta
+
+ betai=1.0/(Rb*remd_t_bath(i))
+ betaiex=1.0/(Rb*remd_t_bath(iex))
+
+ delta=betai*(muca_ene(remd_ene(iex),i,remd_t_bath)-
+ & muca_ene(remd_ene(i),i,remd_t_bath))
+ & -betaiex*(muca_ene(remd_ene(iex),iex,remd_t_bath)-
+ & muca_ene(remd_ene(i),iex,remd_t_bath))
+
+ return
+ end
+
+ double precision function muca_ene(energy,i,remd_t_bath)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ include 'COMMON.MD'
+ double precision y,yp,energy
+ double precision remd_t_bath(maxprocs)
+ integer i
+
+ if (energy.lt.elowi(i)) then
+ call splint(emuca,nemuca,nemuca2,nmuca,elowi(i),y,yp)
+ muca_ene=remd_t_bath(i)*Rb*(yp*(energy-elowi(i))+y)
+ elseif (energy.gt.ehighi(i)) then
+ call splint(emuca,nemuca,nemuca2,nmuca,ehighi(i),y,yp)
+ muca_ene=remd_t_bath(i)*Rb*(yp*(energy-ehighi(i))+y)
+ else
+ call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
+ muca_ene=remd_t_bath(i)*Rb*y
+ endif
+ return
+ end
+
+ subroutine read_muca
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MD'
+ include 'COMMON.REMD'
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
+ imtime=0
+ do i=1,4*maxres
+ hist(i)=0
+ enddo
+ if (modecalc.eq.14.and..not.remd_tlist) then
+ print *,"MUCAREMD works only with TLIST"
+ stop
+ endif
+ open(89,file='muca.input')
+ read(89,*)
+ read(89,*)
+ if (modecalc.eq.14) then
+ read(89,*) (elowi(i),ehighi(i),i=1,nrep)
+ if (remd_mlist) then
+ k=0
+ do i=1,nrep
+ do j=1,remd_m(i)
+ i2rep(k)=i
+ k=k+1
+ enddo
+ enddo
+ elow=elowi(i2rep(me))
+ ehigh=ehighi(i2rep(me))
+ elowi(me+1)=elow
+ ehighi(me+1)=ehigh
+ else
+ elow=elowi(me+1)
+ ehigh=ehighi(me+1)
+ endif
+ else
+ read(89,*) elow,ehigh
+ elowi(1)=elow
+ ehighi(1)=ehigh
+ endif
+ i=0
+ do while(.true.)
+ i=i+1
+ read(89,*,end=100) emuca(i),nemuca(i)
+cd nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
+ enddo
+ 100 continue
+ nmuca=i-1
+ hbin=emuca(nmuca)-emuca(nmuca-1)
+ write (iout,*) 'hbin',hbin
+ write (iout,*) me,'elow,ehigh',elow,ehigh
+ yp1=0
+ ypn=0
+ call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
+ factor_min=0.0d0
+ factor_min=muca_factor(ehigh)
+ call print_muca
+ return
+ end
+
+
+ subroutine print_muca
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MD'
+ include 'COMMON.REMD'
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
+ double precision dummy(maxprocs)
+
+ if (remd_mlist) then
+ k=0
+ do i=1,nrep
+ do j=1,remd_m(i)
+ i2rep(k)=i
+ k=k+1
+ enddo
+ enddo
+ endif
+
+ do i=1,nmuca
+c print *,'nemuca ',emuca(i),nemuca(i)
+ do j=0,4
+ x=emuca(i)+hbin/5*j
+ if (modecalc.eq.14) then
+ if (remd_mlist) then
+ yp=muca_factor(x)*remd_t(i2rep(me))*Rb
+ dummy(me+1)=remd_t(i2rep(me))
+ y=muca_ene(x,me+1,dummy)
+ else
+ yp=muca_factor(x)*remd_t(me+1)*Rb
+ y=muca_ene(x,me+1,remd_t)
+ endif
+ write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
+ & 'muca factor ',x,yp,' muca ene',y
+ else
+ yp=muca_factor(x)*t_bath*Rb
+ dummy(1)=t_bath
+ y=muca_ene(x,1,dummy)
+ write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
+ & 'muca factor ',x,yp,' muca ene',y
+ endif
+ enddo
+ enddo
+ if(mucadyn.gt.0) then
+ do i=1,nmuca
+ write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
+ & imtime,emuca(i),nemuca(i)
+ enddo
+ endif
+ return
+ end
+
+ subroutine muca_update(energy)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MD'
+ include 'COMMON.REMD'
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ double precision energy
+ double precision yp1,ypn
+ integer k
+ logical lnotend
+
+ k=int((energy-emuca(1))/hbin)+1
+
+ IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
+ if(energy.ge.ehigh)
+ & write (iout,*) 'MUCA reject',energy,emuca(k)
+ if(energy.ge.ehigh.and.(energy-ehigh).lt.hbin) then
+ write (iout,*) 'MUCA ehigh',energy,emuca(k)
+ do i=k,nmuca
+ hist(i)=hist(i)+1
+ enddo
+ endif
+ if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1
+ ELSE
+ if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1
+ ENDIF
+ if(mod(imtime,mucadyn).eq.0) then
+
+ do i=1,nmuca
+ IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
+ nemuca(i)=nemuca(i)+dlog(hist(i)+1)
+ ELSE
+ if (hist(i).gt.0) hist(i)=dlog(hist(i))
+ nemuca(i)=nemuca(i)+hist(i)
+ ENDIF
+ hist(i)=0
+ write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
+ & imtime,emuca(i),nemuca(i)
+ enddo
+
+
+ lnotend=.true.
+ ismooth=0
+ ist=2
+ ien=nmuca-1
+ IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
+c lnotend=.false.
+c do i=1,nmuca-1
+c do j=i+1,nmuca
+c if(nemuca(j).lt.nemuca(i)) lnotend=.true.
+c enddo
+c enddo
+ do while(lnotend)
+ ismooth=ismooth+1
+ write (iout,*) 'MUCA update smoothing',ist,ien
+ do i=ist,ien
+ nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
+ enddo
+ lnotend=.false.
+ ist=0
+ ien=0
+ do i=1,nmuca-1
+ do j=i+1,nmuca
+ if(nemuca(j).lt.nemuca(i)) then
+ lnotend=.true.
+ if(ist.eq.0) ist=i-1
+ if(ien.lt.j+1) ien=j+1
+ endif
+ enddo
+ enddo
+ enddo
+ ENDIF
+
+ write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
+ yp1=0
+ ypn=0
+ call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
+ call print_muca
+
+ endif
+ return
+ end
+
+ double precision function muca_factor(energy)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MUCA'
+ double precision y,yp,energy
+
+ if (energy.lt.elow) then
+ call splint(emuca,nemuca,nemuca2,nmuca,elow,y,yp)
+ elseif (energy.gt.ehigh) then
+ call splint(emuca,nemuca,nemuca2,nmuca,ehigh,y,yp)
+ else
+ call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
+ endif
+
+ if(yp.ge.factor_min) then
+ muca_factor=yp
+ else
+ muca_factor=factor_min
+ endif
+cd print *,'energy, muca_factor',energy,muca_factor
+ return
+ end
+
+
+ SUBROUTINE spline(x,y,n,yp1,ypn,y2)
+ INTEGER n,NMAX
+ REAL*8 yp1,ypn,x(n),y(n),y2(n)
+ PARAMETER (NMAX=500)
+ INTEGER i,k
+ REAL*8 p,qn,sig,un,u(NMAX)
+ if (yp1.gt..99e30) then
+ y2(1)=0.
+ u(1)=0.
+ else
+ y2(1)=-0.5
+ u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
+ endif
+ do i=2,n-1
+ sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ p=sig*y2(i-1)+2.
+ y2(i)=(sig-1.)/p
+ u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
+ * /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
+ enddo
+ if (ypn.gt..99e30) then
+ qn=0.
+ un=0.
+ else
+ qn=0.5
+ un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
+ endif
+ y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
+ do k=n-1,1,-1
+ y2(k)=y2(k)*y2(k+1)+u(k)
+ enddo
+ return
+ END
+
+
+ SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
+ INTEGER n
+ REAL*8 x,y,xa(n),y2a(n),ya(n),yp
+ INTEGER k,khi,klo
+ REAL*8 a,b,h
+ klo=1
+ khi=n
+ 1 if (khi-klo.gt.1) then
+ k=(khi+klo)/2
+ if (xa(k).gt.x) then
+ khi=k
+ else
+ klo=k
+ endif
+ goto 1
+ endif
+ h=xa(khi)-xa(klo)
+ if (h.eq.0.) pause 'bad xa input in splint'
+ a=(xa(khi)-x)/h
+ b=(x-xa(klo))/h
+ y=a*ya(klo)+b*ya(khi)+
+ * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
+ yp=-ya(klo)/h+ya(khi)/h-3*(a**2)*y2a(klo)*h/6.
+ + +(3*(b**2)-1)*y2a(khi)*h/6.
+ return
+ END
projects / unres.git / blobdiff