+++ /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