主题:[讨论]求救 !!!!!!!!!!利用 蒙地卡罗计算 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!請問大大哪裡錯了
计算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!請問大大哪裡錯了