回 帖 发 新 帖 刷新版面

主题:[讨论]求救 !!!!!!!!!!利用 蒙地卡罗计算 butane

利用 蒙地卡罗
计算GAUCHE FORM(角度大于等于50小于等于70度)
ANTI FORM(角度角度大于等于170小于等于190度)
不同比例下 比例与温度之关系!
T为角度
E(T)=v1+v1*cost+v2-v2*cos2t
+v3+v3*cos3t

v1=1.25,v2=0.025   v3=1.5

(p=exp(-deltaE/RT)
R=1.9858775E-3
我不会写.................可以请各位高手一步一步教我吗~~~QQ…谢谢……………….


program main
real R,THERMALE,TEMP,angle,ang_deg,energy,engx,Nanti,Ngau,pi
integer nsearch
write(*,*)'input the number of times to search!!!'
read(*,*)nsearch
write(*,*)'input the temperature in K!!!'
read(*,*)temp
pi=3.14159 engx=10
nanti=0
Ngau=0
R=1.9858775E-3
THERMALE=R*temp
do I=1,36
confomers(I)=0
enddo
write(*,*)'degree energy ran# prob,engx'
do I=1,Nsearch

x=ran( )
ang_deg=x*360.0
angle=(ang_deg)*pi/180.0
call Butan(I,angle,energy)
If (energy.lt.engx)then
write(3,100)'accepted!',ang_deg,energy,x
engx=energy
N=NINT(ang_deg/10.0)
confomers(N)=confomers(N)+1 if(ang_deg.ge.170.0)and(ang_deg.le.190.0)Nanti=nanti+1

if(ang_deg.ge.50.0)and(ang_deg.le.70.0)Ngau=Ngau+1
else
x=ran()
prob=exp(-(energy-engx)/THERMALE)
If (prob.ge.x)then
write(3,100)'Probility!',ang_deg,energy,x,prob
engx=Energy
N=NINT(ang_deg/10.0)
conformers(N)=conformers(N)+1
else

write(3,120)'rejected!',ang_deg,energy,x,prob,engx
Endif
endif
enddo

do I=1,36
write(*,*)'Torision Angle ',(I-1)*10,'-',I*10,conformers(I)
enddo
write(*,130)'# of anti-forms=',Nanti
write(*,130)'# of gauche-forms=',ngau
write(*,130)'# of Gauche/anti ratio=',Ngau/Nanti
100format(a11,3,(2x,f10.5))

120format(a11,4,(2x,f10.5))
130format(a11,5,(2x,f10.5))

131 format(a25,1X,F8.4)

140 format(a25,1x,I4,A3,I4,1x,i6)


close(3)
stop
end
subroutine Butan(J,ang,eng)
integer J
real ang,ang2,ang3,eng,v0,v1,v2,v3
ang2=ang*2.0
ang3=ang*3.0
Eng=v0+v1+v2+v3+v1*cos(ang)-V2*cos(ang2)+v3*cos(ang3)
200 format(i6,2x,F8.5,2X,F8.5,2X,F8.5)
return
end
以上很多要DEBUG!請問大大哪裡錯了

回复列表 (共2个回复)

沙发

请各位老师前辈帮帮忙啊~救救小妹!!

板凳

蒙特卡洛或者蒙特卡罗吧

我来回复

您尚未登录,请登录后再回复。点此登录或注册