module parameters implicit none integer, parameter :: long = selected_real_kind(15,307) integer, parameter :: na=1680,nrep=24 real(long), parameter :: RT=0.596d0, temp0=300.d0 real(long), dimension(na,3) :: vel1,vel2 real(long), dimension(nrep) :: temp,energy integer:: mseed character(len=40) :: file1,file2 end module parameters program exchange use parameters implicit none integer :: i,m,n,run,repstep,l real(long) :: ran2,scale,r,expd,delta logical :: Accept read*,repstep,run mseed=-run*repstep - 51000 open(20,file='RepTemp.dat') do i=1,nrep read(20,*)temp(i) temp(i) = RT*temp(i)/temp0 enddo close(20) open(20,file='RepEnergy.dat') do i=1,nrep read(20,*)energy(i) enddo close(20) l = int(2*ran2(mseed)) + 1 do m=l,nrep-1,2 Accept=.false. n = m+1 delta = (energy(m)-energy(n))*(1.d0/temp(n)-1.d0/temp(m)) if(delta <= 0.d0) then Accept=.true. else r = ran2(mseed) expd = exp(-delta) if(r < expd) Accept=.true. endif if(Accept)then print*,m,n,1 if(m < 10) file1 = 'velocities'//char(48+m) if(m >= 10) file1 = 'velocities'//char(48+m/10)//char(48+m-m/10*10) if(n < 10) file2 = 'velocities'//char(48+n) if(n >= 10) file2 = 'velocities'//char(48+n/10)//char(48+n-n/10*10) open(21,file=trim(file1)//'.dat') open(22,file=trim(file2)//'.dat') do i=1,na read(21,100)vel1(i,:) read(22,100)vel2(i,:) enddo close(21) close(22) scale=sqrt(temp(n)/temp(m)) open(21,file=trim(file1)//'s.dat') open(22,file=trim(file2)//'s.dat') do i=1,na write(21,100)vel2(i,:)/scale write(22,100)scale*vel1(i,:) enddo close(21) close(22) else print*,m,n,0 endif enddo 100 format(3d22.15) end program exchange !-------------------------------------------------------------------- function ran2(idum) use parameters, only : long implicit none integer, intent(inout) :: idum integer, parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1, & IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, & NTAB=32,NDIV=1+IMM1/NTAB real (long), parameter :: AM=1.D0/IM1,ESP1=1.2D-7,RNMX=1.D0-ESP1 integer :: idum2,j,k,iy integer, dimension(NTAB) :: iv real (long) :: ran2 save iv,iy,idum2 data idum2/123456789/, iv/NTAB*0/, iy/0/ if(idum <= 0) then idum=max(-idum,1) idum2=idum do j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum < 0) idum=idum+IM1 if (j <= NTAB) iv(j)=idum enddo iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum < 0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if(idum2 < 0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy < 1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) end function ran2