1 subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
8 include 'COMMON.LANGEVIN.lang0.5diag'
10 include 'COMMON.LANGEVIN.lang0'
13 include 'COMMON.LANGEVIN'
15 double precision remd_t_bath(maxprocs)
16 double precision remd_ene(maxprocs)
17 double precision muca_ene
18 double precision betai,betaiex,delta
21 betai=1.0/(Rb*remd_t_bath(i))
22 betaiex=1.0/(Rb*remd_t_bath(iex))
24 delta=betai*(muca_ene(remd_ene(iex),i,remd_t_bath)-
25 & muca_ene(remd_ene(i),i,remd_t_bath))
26 & -betaiex*(muca_ene(remd_ene(iex),iex,remd_t_bath)-
27 & muca_ene(remd_ene(i),iex,remd_t_bath))
32 double precision function muca_ene(energy,i,remd_t_bath)
39 include 'COMMON.LANGEVIN.lang0.5diag'
41 include 'COMMON.LANGEVIN.lang0'
44 include 'COMMON.LANGEVIN'
46 double precision y,yp,energy
47 double precision remd_t_bath(maxprocs)
50 if (energy.lt.elowi(i)) then
51 call splint(emuca,nemuca,nemuca2,nmuca,elowi(i),y,yp)
52 muca_ene=remd_t_bath(i)*Rb*(yp*(energy-elowi(i))+y)
53 elseif (energy.gt.ehighi(i)) then
54 call splint(emuca,nemuca,nemuca2,nmuca,ehighi(i),y,yp)
55 muca_ene=remd_t_bath(i)*Rb*(yp*(energy-ehighi(i))+y)
57 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
58 muca_ene=remd_t_bath(i)*Rb*y
67 include 'COMMON.CONTROL'
70 include 'COMMON.SETUP'
71 include 'COMMON.IOUNITS'
72 double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
78 if (modecalc.eq.14.and..not.remd_tlist) then
79 print *,"MUCAREMD works only with TLIST"
82 open(89,file='muca.input')
85 if (modecalc.eq.14) then
86 read(89,*) (elowi(i),ehighi(i),i=1,nrep)
96 ehigh=ehighi(i2rep(me))
104 read(89,*) elow,ehigh
111 read(89,*,end=100) emuca(i),nemuca(i)
112 cd nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
116 hbin=emuca(nmuca)-emuca(nmuca-1)
117 write (iout,*) 'hbin',hbin
118 write (iout,*) me,'elow,ehigh',elow,ehigh
121 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
123 factor_min=muca_factor(ehigh)
129 subroutine print_muca
132 include 'COMMON.MUCA'
133 include 'COMMON.CONTROL'
137 include 'COMMON.LANGEVIN.lang0.5diag'
139 include 'COMMON.LANGEVIN.lang0'
142 include 'COMMON.LANGEVIN'
144 include 'COMMON.REMD'
145 include 'COMMON.SETUP'
146 include 'COMMON.IOUNITS'
147 double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
148 double precision dummy(maxprocs)
161 c print *,'nemuca ',emuca(i),nemuca(i)
164 if (modecalc.eq.14) then
166 yp=muca_factor(x)*remd_t(i2rep(me))*Rb
167 dummy(me+1)=remd_t(i2rep(me))
168 y=muca_ene(x,me+1,dummy)
170 yp=muca_factor(x)*remd_t(me+1)*Rb
171 y=muca_ene(x,me+1,remd_t)
173 write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
174 & 'muca factor ',x,yp,' muca ene',y
176 yp=muca_factor(x)*t_bath*Rb
178 y=muca_ene(x,1,dummy)
179 write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
180 & 'muca factor ',x,yp,' muca ene',y
184 if(mucadyn.gt.0) then
186 write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
187 & imtime,emuca(i),nemuca(i)
193 subroutine muca_update(energy)
196 include 'COMMON.MUCA'
197 include 'COMMON.CONTROL'
199 include 'COMMON.REMD'
200 include 'COMMON.SETUP'
201 include 'COMMON.IOUNITS'
202 double precision energy
203 double precision yp1,ypn
204 integer i,j,k,ismooth,ist,ien
207 k=int((energy-emuca(1))/hbin)+1
209 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
211 & write (iout,*) 'MUCA reject',energy,emuca(k)
212 if(energy.ge.ehigh.and.(energy-ehigh).lt.hbin) then
213 write (iout,*) 'MUCA ehigh',energy,emuca(k)
218 if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1
220 if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1
222 if(mod(imtime,mucadyn).eq.0) then
225 IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
226 nemuca(i)=nemuca(i)+dlog(hist(i)+1)
228 if (hist(i).gt.0) hist(i)=dlog(hist(i))
229 nemuca(i)=nemuca(i)+hist(i)
232 write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
233 & imtime,emuca(i),nemuca(i)
241 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
245 c if(nemuca(j).lt.nemuca(i)) lnotend=.true.
250 write (iout,*) 'MUCA update smoothing',ist,ien
252 nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
259 if(nemuca(j).lt.nemuca(i)) then
262 if(ien.lt.j+1) ien=j+1
269 write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
272 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
279 double precision function muca_factor(energy)
282 include 'COMMON.MUCA'
283 double precision y,yp,energy
285 if (energy.lt.elow) then
286 call splint(emuca,nemuca,nemuca2,nmuca,elow,y,yp)
287 elseif (energy.gt.ehigh) then
288 call splint(emuca,nemuca,nemuca2,nmuca,ehigh,y,yp)
290 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
293 if(yp.ge.factor_min) then
296 muca_factor=factor_min
298 cd print *,'energy, muca_factor',energy,muca_factor
303 SUBROUTINE spline(x,y,n,yp1,ypn,y2)
306 REAL*8 yp1,ypn,x(n),y(n),y2(n)
309 REAL*8 p,qn,sig,un,u(NMAX)
310 if (yp1.gt..99e30) then
315 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
318 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
321 u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
322 * /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
324 if (ypn.gt..99e30) then
329 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
331 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
333 y2(k)=y2(k)*y2(k+1)+u(k)
339 SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
342 REAL*8 x,y,xa(n),y2a(n),ya(n),yp
347 1 if (khi-klo.gt.1) then
357 if (h.eq.0.) pause 'bad xa input in splint'
360 y=a*ya(klo)+b*ya(khi)+
361 * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
362 yp=-ya(klo)/h+ya(khi)/h-3*(a**2)*y2a(klo)*h/6.
363 + +(3*(b**2)-1)*y2a(khi)*h/6.