1 subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
2 implicit real*8 (a-h,o-z)
6 double precision remd_t_bath(maxprocs)
7 double precision remd_ene(maxprocs)
8 double precision muca_ene
9 double precision betai,betaiex,delta
11 betai=1.0/(Rb*remd_t_bath(i))
12 betaiex=1.0/(Rb*remd_t_bath(iex))
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))
22 double precision function muca_ene(energy,i,remd_t_bath)
23 implicit real*8 (a-h,o-z)
27 double precision y,yp,energy
28 double precision remd_t_bath(maxprocs)
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)
38 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
39 muca_ene=remd_t_bath(i)*Rb*y
45 implicit real*8 (a-h,o-z)
48 include 'COMMON.CONTROL'
51 include 'COMMON.SETUP'
52 include 'COMMON.IOUNITS'
53 double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
58 if (modecalc.eq.14.and..not.remd_tlist) then
59 print *,"MUCAREMD works only with TLIST"
62 open(89,file='muca.input')
65 if (modecalc.eq.14) then
66 read(89,*) (elowi(i),ehighi(i),i=1,nrep)
76 ehigh=ehighi(i2rep(me))
91 read(89,*,end=100) emuca(i),nemuca(i)
92 cd nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
96 hbin=emuca(nmuca)-emuca(nmuca-1)
97 write (iout,*) 'hbin',hbin
98 write (iout,*) me,'elow,ehigh',elow,ehigh
101 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
103 factor_min=muca_factor(ehigh)
109 subroutine print_muca
110 implicit real*8 (a-h,o-z)
112 include 'COMMON.MUCA'
113 include 'COMMON.CONTROL'
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)
132 c print *,'nemuca ',emuca(i),nemuca(i)
135 if (modecalc.eq.14) 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)
141 yp=muca_factor(x)*remd_t(me+1)*Rb
142 y=muca_ene(x,me+1,remd_t)
144 write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
145 & 'muca factor ',x,yp,' muca ene',y
147 yp=muca_factor(x)*t_bath*Rb
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
155 if(mucadyn.gt.0) then
157 write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
158 & imtime,emuca(i),nemuca(i)
164 subroutine muca_update(energy)
165 implicit real*8 (a-h,o-z)
167 include 'COMMON.MUCA'
168 include 'COMMON.CONTROL'
170 include 'COMMON.REMD'
171 include 'COMMON.SETUP'
172 include 'COMMON.IOUNITS'
173 double precision energy
174 double precision yp1,ypn
178 k=int((energy-emuca(1))/hbin)+1
180 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
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)
189 if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1
191 if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1
193 if(mod(imtime,mucadyn).eq.0) then
196 IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
197 nemuca(i)=nemuca(i)+dlog(hist(i)+1)
199 if (hist(i).gt.0) hist(i)=dlog(hist(i))
200 nemuca(i)=nemuca(i)+hist(i)
203 write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
204 & imtime,emuca(i),nemuca(i)
212 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
216 c if(nemuca(j).lt.nemuca(i)) lnotend=.true.
221 write (iout,*) 'MUCA update smoothing',ist,ien
223 nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
230 if(nemuca(j).lt.nemuca(i)) then
233 if(ien.lt.j+1) ien=j+1
240 write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
243 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
250 double precision function muca_factor(energy)
251 implicit real*8 (a-h,o-z)
253 include 'COMMON.MUCA'
254 double precision y,yp,energy
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)
261 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
264 if(yp.ge.factor_min) then
267 muca_factor=factor_min
269 cd print *,'energy, muca_factor',energy,muca_factor
274 SUBROUTINE spline(x,y,n,yp1,ypn,y2)
276 REAL*8 yp1,ypn,x(n),y(n),y2(n)
279 REAL*8 p,qn,sig,un,u(NMAX)
280 if (yp1.gt..99e30) then
285 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
288 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
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
294 if (ypn.gt..99e30) then
299 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
301 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
303 y2(k)=y2(k)*y2(k+1)+u(k)
309 SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
311 REAL*8 x,y,xa(n),y2a(n),ya(n),yp
316 1 if (khi-klo.gt.1) then
326 if (h.eq.0.) pause 'bad xa input in splint'
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.