2 !-----------------------------------------------------------------------------
6 use geometry_data, only:nres
8 !-----------------------------------------------------------------------------
10 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
13 subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
16 ! implicit real*8 (a-h,o-z)
17 ! include 'DIMENSIONS'
18 ! include 'COMMON.MUCA'
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
26 betai=1.0/(Rb*remd_t_bath(i))
27 betaiex=1.0/(Rb*remd_t_bath(iex))
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))
35 end subroutine muca_delta
36 !-----------------------------------------------------------------------------
37 real(kind=8) function muca_ene(energy,i,remd_t_bath)
40 ! implicit real*8 (a-h,o-z)
41 ! include 'DIMENSIONS'
42 ! include 'COMMON.MUCA'
44 real(kind=8) :: y,yp,energy
45 real(kind=8),dimension(nprocs) :: remd_t_bath !(maxprocs)
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)
55 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
56 muca_ene=remd_t_bath(i)*Rb*y
60 !-----------------------------------------------------------------------------
61 subroutine muca_update(energy)
64 ! implicit real*8 (a-h,o-z)
65 ! include 'DIMENSIONS'
66 ! include 'COMMON.MUCA'
67 ! include 'COMMON.CONTROL'
69 ! include 'COMMON.REMD'
70 ! include 'COMMON.SETUP'
71 ! include 'COMMON.IOUNITS'
72 real(kind=8) :: energy
73 real(kind=8) :: yp1,ypn
74 integer :: k,i,ismooth,ist,ien,j
77 k=int((energy-emuca(1))/hbin)+1
79 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
81 write (iout,*) 'MUCA reject',energy,emuca(k)
82 if(energy.ge.ehigh.and.(energy-ehigh).lt.hbin) then
83 write (iout,*) 'MUCA ehigh',energy,emuca(k)
88 if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1
90 if(k.gt.0.and.k.lt.4*nres) hist(k)=hist(k)+1
92 if(mod(imtime,mucadyn).eq.0) then
95 IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
96 nemuca(i)=nemuca(i)+dlog(hist(i)+1)
98 if (hist(i).gt.0) hist(i)=dlog(hist(i))
99 nemuca(i)=nemuca(i)+hist(i)
102 write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',&
103 imtime,emuca(i),nemuca(i)
111 IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
115 ! if(nemuca(j).lt.nemuca(i)) lnotend=.true.
120 write (iout,*) 'MUCA update smoothing',ist,ien
122 nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
129 if(nemuca(j).lt.nemuca(i)) then
132 if(ien.lt.j+1) ien=j+1
139 write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
142 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
147 end subroutine muca_update
148 !-----------------------------------------------------------------------------
149 real(kind=8) function muca_factor(energy)
151 ! implicit real*8 (a-h,o-z)
152 ! include 'DIMENSIONS'
153 ! include 'COMMON.MUCA'
154 real(kind=8) :: y,yp,energy
156 if (energy.lt.elow) then
157 call splint(emuca,nemuca,nemuca2,nmuca,elow,y,yp)
158 elseif (energy.gt.ehigh) then
159 call splint(emuca,nemuca,nemuca2,nmuca,ehigh,y,yp)
161 call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
164 if(yp.ge.factor_min) then
167 muca_factor=factor_min
169 !d print *,'energy, muca_factor',energy,muca_factor
171 end function muca_factor
172 !-----------------------------------------------------------------------------
173 subroutine spline(x,y,n,yp1,ypn,y2)
176 REAL(kind=8) :: yp1,ypn,x(n),y(n),y2(n)
177 integer,PARAMETER :: NMAX=500
179 REAL(kind=8) :: p,qn,sig,un,u(NMAX)
180 if (yp1.gt..99e30) then
185 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
188 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
191 u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
192 /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
194 if (ypn.gt..99e30) then
199 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
201 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
203 y2(k)=y2(k)*y2(k+1)+u(k)
206 end subroutine spline
207 !-----------------------------------------------------------------------------
208 subroutine splint(xa,ya,y2a,n,x,y,yp)
211 REAL(kind=8) :: x,y,xa(n),y2a(n),ya(n),yp
213 REAL(kind=8) :: a,b,h
216 1 if (khi-klo.gt.1) then
226 if (h.eq.0.) pause 'bad xa input in splint'
229 y=a*ya(klo)+b*ya(khi)+ &
230 ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
231 yp=-ya(klo)/h+ya(khi)/h-3*(a**2)*y2a(klo)*h/6. &
232 +(3*(b**2)-1)*y2a(khi)*h/6.
234 end subroutine splint
235 !-----------------------------------------------------------------------------
237 !-----------------------------------------------------------------------------
240 use control_data, only: modecalc !,maxprocs
244 ! implicit real*8 (a-h,o-z)
245 ! include 'DIMENSIONS'
246 ! include 'COMMON.MUCA'
247 ! include 'COMMON.CONTROL'
248 ! include 'COMMON.MD'
249 ! include 'COMMON.REMD'
250 ! include 'COMMON.SETUP'
251 ! include 'COMMON.IOUNITS'
252 real(kind=8) :: yp1,ypn,yp,x,y,muca_ene !,muca_factor
254 real(kind=8) :: emuca_alloc(4*maxres),nemuca_alloc(4*maxres)
255 integer(kind=2) :: i2rep_alloc(0:Nprocs) !(0:maxprocs)
257 ! real(kind=8) :: var,ene
260 allocate(hist(4*maxres)) !(4*maxres)
261 !el allocate(i2rep(0:nodes+1)) !(0:maxprocs)
265 if (modecalc.eq.14.and..not.remd_tlist) then
266 print *,"MUCAREMD works only with TLIST"
269 open(89,file='muca.input')
272 if (modecalc.eq.14) then
273 allocate(elowi(nrep),ehighi(nrep)) !(maxprocs)
274 read(89,*) (elowi(i),ehighi(i),i=1,nrep)
284 allocate(i2rep(k)) !(0:maxprocs)
286 i2rep(j)=i2rep_alloc(j)
288 elow=elowi(i2rep(me))
289 ehigh=ehighi(i2rep(me))
297 read(89,*) elow,ehigh
304 read(89,*,end=100) emuca_alloc(i),nemuca_alloc(i)
305 !d nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
307 allocate(emuca(i),nemuca(i),nemuca2(i)) !4*maxres
309 emuca(j)=emuca_alloc(j)
310 nemuca(j)=nemuca_alloc(j)
314 hbin=emuca(nmuca)-emuca(nmuca-1)
315 write (iout,*) 'hbin',hbin
316 write (iout,*) me,'elow,ehigh',elow,ehigh
319 call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
321 factor_min=muca_factor(ehigh)
324 end subroutine read_muca
325 !-----------------------------------------------------------------------------
326 subroutine print_muca
328 use control_data, only: modecalc,mucadyn !,maxprocs
332 ! implicit real*8 (a-h,o-z)
333 ! include 'DIMENSIONS'
334 ! include 'COMMON.MUCA'
335 ! include 'COMMON.CONTROL'
336 ! include 'COMMON.MD'
337 ! include 'COMMON.REMD'
338 ! include 'COMMON.SETUP'
339 ! include 'COMMON.IOUNITS'
340 real(kind=8) :: yp1,ypn,yp,x,y !,muca_ene,muca_factor
341 real(kind=8) :: dummy(Nprocs) !(maxprocs)
356 ! print *,'nemuca ',emuca(i),nemuca(i)
359 if (modecalc.eq.14) then
361 yp=muca_factor(x)*remd_t(i2rep(me))*Rb
362 dummy(me+1)=remd_t(i2rep(me))
363 y=muca_ene(x,me+1,dummy)
365 yp=muca_factor(x)*remd_t(me+1)*Rb
366 y=muca_ene(x,me+1,remd_t)
368 write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,&
369 'muca factor ',x,yp,' muca ene',y
371 yp=muca_factor(x)*t_bath*Rb
373 y=muca_ene(x,1,dummy)
374 write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,&
375 'muca factor ',x,yp,' muca ene',y
379 if(mucadyn.gt.0) then
381 write(iout,'(a13,i8,2f12.5)') 'nemuca after ',&
382 imtime,emuca(i),nemuca(i)
386 end subroutine print_muca
387 !-----------------------------------------------------------------------------
388 !-----------------------------------------------------------------------------