X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fsrc%2Fmuca_md.f;fp=source%2Funres%2Fsrc_MD%2Fsrc%2Fmuca_md.f;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=c10a6a763c86869fc72689d489252476d83c4198;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/src/muca_md.f b/source/unres/src_MD/src/muca_md.f deleted file mode 100644 index c10a6a7..0000000 --- a/source/unres/src_MD/src/muca_md.f +++ /dev/null @@ -1,334 +0,0 @@ - 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