c10a6a763c86869fc72689d489252476d83c4198
[unres.git] / source / unres / src_MD / src / muca_md.f
1       subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
2       implicit real*8 (a-h,o-z)     
3       include 'DIMENSIONS'
4       include 'COMMON.MUCA'
5       include 'COMMON.MD'
6       double precision remd_t_bath(maxprocs)
7       double precision remd_ene(maxprocs)
8       double precision muca_ene
9       double precision betai,betaiex,delta
10
11       betai=1.0/(Rb*remd_t_bath(i))
12       betaiex=1.0/(Rb*remd_t_bath(iex))
13       
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))
18
19       return
20       end
21       
22       double precision function muca_ene(energy,i,remd_t_bath)
23       implicit real*8 (a-h,o-z)
24       include 'DIMENSIONS'
25       include 'COMMON.MUCA'
26       include 'COMMON.MD'
27       double precision y,yp,energy
28       double precision remd_t_bath(maxprocs)
29       integer i
30  
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)
37       else
38         call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
39         muca_ene=remd_t_bath(i)*Rb*y
40       endif
41       return
42       end
43       
44       subroutine read_muca
45       implicit real*8 (a-h,o-z)
46       include 'DIMENSIONS'
47       include 'COMMON.MUCA'
48       include 'COMMON.CONTROL'
49       include 'COMMON.MD'
50       include 'COMMON.REMD'
51       include 'COMMON.SETUP'
52       include 'COMMON.IOUNITS'
53       double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
54       imtime=0
55       do i=1,4*maxres
56         hist(i)=0
57       enddo
58       if (modecalc.eq.14.and..not.remd_tlist) then
59                 print *,"MUCAREMD works only with TLIST"
60                 stop
61       endif
62       open(89,file='muca.input')      
63       read(89,*) 
64       read(89,*) 
65       if (modecalc.eq.14) then
66         read(89,*) (elowi(i),ehighi(i),i=1,nrep)
67        if (remd_mlist) then
68         k=0
69         do i=1,nrep
70          do j=1,remd_m(i)
71           i2rep(k)=i
72           k=k+1
73          enddo
74         enddo
75         elow=elowi(i2rep(me))
76         ehigh=ehighi(i2rep(me))
77         elowi(me+1)=elow
78         ehighi(me+1)=ehigh        
79        else 
80         elow=elowi(me+1)
81         ehigh=ehighi(me+1)
82        endif
83       else
84         read(89,*) elow,ehigh
85         elowi(1)=elow
86         ehighi(1)=ehigh
87       endif
88       i=0
89       do while(.true.)
90        i=i+1
91        read(89,*,end=100) emuca(i),nemuca(i)
92 cd      nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
93       enddo
94  100  continue
95       nmuca=i-1
96       hbin=emuca(nmuca)-emuca(nmuca-1)
97       write (iout,*) 'hbin',hbin      
98       write (iout,*) me,'elow,ehigh',elow,ehigh
99       yp1=0
100       ypn=0
101       call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
102       factor_min=0.0d0
103       factor_min=muca_factor(ehigh)
104       call print_muca
105       return
106       end
107
108
109       subroutine print_muca
110       implicit real*8 (a-h,o-z)
111       include 'DIMENSIONS'
112       include 'COMMON.MUCA'
113       include 'COMMON.CONTROL'
114       include 'COMMON.MD'
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)
120
121       if (remd_mlist) then
122            k=0
123            do i=1,nrep
124             do j=1,remd_m(i)
125              i2rep(k)=i
126              k=k+1
127             enddo
128            enddo
129       endif
130
131       do i=1,nmuca
132 c       print *,'nemuca ',emuca(i),nemuca(i)
133        do j=0,4
134         x=emuca(i)+hbin/5*j
135         if (modecalc.eq.14) then
136          if (remd_mlist) 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)
140          else
141           yp=muca_factor(x)*remd_t(me+1)*Rb
142           y=muca_ene(x,me+1,remd_t)
143          endif
144          write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
145      &      'muca factor ',x,yp,' muca ene',y
146         else
147          yp=muca_factor(x)*t_bath*Rb
148          dummy(1)=t_bath
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         
152         endif
153        enddo
154       enddo
155       if(mucadyn.gt.0) then
156        do i=1,nmuca
157          write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
158      &             imtime,emuca(i),nemuca(i)
159        enddo
160       endif
161       return
162       end
163
164       subroutine muca_update(energy)
165       implicit real*8 (a-h,o-z)
166       include 'DIMENSIONS'
167       include 'COMMON.MUCA'
168       include 'COMMON.CONTROL'
169       include 'COMMON.MD'
170       include 'COMMON.REMD'
171       include 'COMMON.SETUP'
172       include 'COMMON.IOUNITS'
173       double precision energy
174       double precision yp1,ypn
175       integer k
176       logical lnotend
177
178       k=int((energy-emuca(1))/hbin)+1
179       
180       IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
181        if(energy.ge.ehigh) 
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)
185          do i=k,nmuca
186            hist(i)=hist(i)+1
187          enddo
188        endif
189        if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1       
190       ELSE
191        if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1       
192       ENDIF
193       if(mod(imtime,mucadyn).eq.0) then
194
195          do i=1,nmuca
196           IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
197            nemuca(i)=nemuca(i)+dlog(hist(i)+1)
198           ELSE
199            if (hist(i).gt.0) hist(i)=dlog(hist(i))         
200            nemuca(i)=nemuca(i)+hist(i)
201           ENDIF
202           hist(i)=0
203           write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
204      &          imtime,emuca(i),nemuca(i)
205          enddo
206
207
208          lnotend=.true.
209          ismooth=0
210          ist=2
211          ien=nmuca-1
212         IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN         
213 c         lnotend=.false.         
214 c         do i=1,nmuca-1
215 c           do j=i+1,nmuca
216 c            if(nemuca(j).lt.nemuca(i)) lnotend=.true.
217 c           enddo
218 c         enddo         
219          do while(lnotend)
220           ismooth=ismooth+1
221           write (iout,*) 'MUCA update smoothing',ist,ien
222           do i=ist,ien
223            nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
224           enddo
225           lnotend=.false.
226           ist=0
227           ien=0
228           do i=1,nmuca-1
229            do j=i+1,nmuca
230             if(nemuca(j).lt.nemuca(i)) then 
231               lnotend=.true.
232               if(ist.eq.0) ist=i-1
233               if(ien.lt.j+1) ien=j+1
234             endif
235            enddo
236           enddo
237          enddo
238         ENDIF 
239
240          write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
241          yp1=0
242          ypn=0
243          call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
244          call print_muca
245          
246       endif
247       return
248       end
249       
250       double precision function muca_factor(energy)
251       implicit real*8 (a-h,o-z)
252       include 'DIMENSIONS'
253       include 'COMMON.MUCA'
254       double precision y,yp,energy
255       
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)
260       else
261         call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
262       endif
263       
264       if(yp.ge.factor_min) then
265        muca_factor=yp
266       else
267        muca_factor=factor_min
268       endif
269 cd      print *,'energy, muca_factor',energy,muca_factor
270       return
271       end
272       
273
274       SUBROUTINE spline(x,y,n,yp1,ypn,y2)
275       INTEGER n,NMAX
276       REAL*8 yp1,ypn,x(n),y(n),y2(n)
277       PARAMETER (NMAX=500)
278       INTEGER i,k
279       REAL*8 p,qn,sig,un,u(NMAX)
280       if (yp1.gt..99e30) then 
281       y2(1)=0. 
282       u(1)=0.
283       else 
284          y2(1)=-0.5
285          u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
286       endif
287       do i=2,n-1 
288          sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
289          p=sig*y2(i-1)+2.
290          y2(i)=(sig-1.)/p
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
293       enddo 
294       if (ypn.gt..99e30) then 
295          qn=0.
296          un=0.
297       else 
298          qn=0.5
299          un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
300       endif
301       y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
302       do k=n-1,1,-1 
303          y2(k)=y2(k)*y2(k+1)+u(k) 
304       enddo 
305       return
306       END 
307
308
309       SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
310       INTEGER n
311       REAL*8 x,y,xa(n),y2a(n),ya(n),yp
312       INTEGER k,khi,klo
313       REAL*8 a,b,h
314       klo=1 
315       khi=n
316  1    if (khi-klo.gt.1) then
317          k=(khi+klo)/2
318          if (xa(k).gt.x) then
319             khi=k
320          else
321             klo=k
322          endif
323          goto 1
324       endif 
325       h=xa(khi)-xa(klo)
326       if (h.eq.0.) pause 'bad xa input in splint' 
327       a=(xa(khi)-x)/h 
328       b=(x-xa(klo))/h
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.
333       return
334       END