X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2FSRC-SURPLUS%2Fmuca_md.f;fp=source%2Funres%2Fsrc_Eshel%2FSRC-SURPLUS%2Fmuca_md.f;h=c10a6a763c86869fc72689d489252476d83c4198;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_Eshel/SRC-SURPLUS/muca_md.f b/source/unres/src_Eshel/SRC-SURPLUS/muca_md.f new file mode 100644 index 0000000..c10a6a7 --- /dev/null +++ b/source/unres/src_Eshel/SRC-SURPLUS/muca_md.f @@ -0,0 +1,334 @@ + 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