| 查看: 3174 | 回复: 7 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
Gina88木虫 (正式写手)
|
[求助]
Brich-Murnaghan状态方程去拟合平衡时的晶格体积出问题！
|
|
|
用Brich-Murnaghan状态方程去拟合平衡时的晶格体积。用体积回代后的值怎么和真实值差那么多啊？ !announcement begin******************** implicit none real(8),allocatable::V(  real(8),allocatable::E( real(8),allocatable::EV( real(8),allocatable::Ecal( real(8)::b(4),x(4) real(8)::a(4,4) real(8)::sum real(8)::E0,V0,B0,B01 real(8)::V0_2,B0_2,B01_2 integer i,j,k,n,ios character(len=80)::nonsense !End announcement********************** !Begin reading data from "comment_EvsV"******** open(unit=18,file="comment_EvsV",status="old"  n=0 do read(18,*,iostat=ios)nonsense if (ios < 0) exit n=n+1 end do close(18) allocate(V(n)) allocate(E(n)) allocate(EV(n)) open(unit=18,file="comment_EvsV",status="old" do i=1,n read(18,*)V(i),E(i) end do close(18) !End reading data from "comment_EvsV"********** do i=1,4 do j=1,4 b(i)=0 a(i,j)=0 end do end do do i=1,n a(1,1)=a(1,1)+V(i)**(-4) a(1,2)=a(1,2)+V(i)**(-10/3) a(1,3)=a(1,3)+V(i)**(-8/3) a(1,4)=a(1,4)+V(i)**(-2) a(2,4)=a(2,4)+V(i)**(-4/3) a(3,4)=a(3,4)+V(i)**(-2/3) b(1)=b(1)+E(i)*(V(i)**(-2)) b(2)=b(2)+E(i)*(V(i)**(-4/3)) b(3)=b(3)+E(i)*(V(i)**(-2/3)) b(4)=b(4)+E(i) end do a(2,1)=a(1,2) a(2,2)=a(1,3);a(3,1)=a(1,3) a(2,3)=a(1,4);a(3,2)=a(1,4);a(4,1)=a(1,4) a(3,3)=a(2,4);a(4,2)=a(2,4) a(4,3)=a(3,4) a(4,4)=n do k=1,3 !if(a(k,k)=0)then ! print(*,'(a)')"divide by 0" ! exit !end if do i=k+1,4 do j=k+1,4 a(i,j)=a(i,j)+a(k,j)*(-a(i,k)/a(k,k)) b(i)=b(i)+b(k)*(-a(i,k)/a(k,k)) end do end do end do x(4)=b(4)/a(4,4) do i=3,1,-1 sum=0 do j=i+1,4 sum=sum+a(i,j)*x(j) end do x(i)=(b(i)-sum)/a(i,i) end do E0=x(4) sum=(42*x(1)*x(3)/(x(2)**2)-14)**2 sum=sum-(9*x(1)*x(3)/(x(2)**2)-3)*(196*x(1)*x(3)/(x(2)**2)-64) B01=((42*x(1)*x(3)/(x(2)**2)-14)+sqrt(sum))/(9*x(1)*x(3)/(x(2)**2)-3) B01_2=((42*x(1)*x(3)/(x(2)**2)-14)-sqrt(sum))/(9*x(1)*x(3)/(x(2)**2)-3) sum=-(x(1)/x(2))*((3*B01-14)/(B01-4)) if(sum>0)then V0=sum**1.5 B0=(16/9)*x(2)*V0**(-7/3)/(14-3*B01) sum=-(x(1)/x(2))*((3*B01_2-14)/(B01_2-4)) else if (sum>0)then V0_2=sum**1.5 B0_2=(16/9)*x(2)*V0_2**(-7/3)/(14-3*B01_2) end if allocate(Ecal(n)) do i=1,n Ecal(i)=x(1)*(V(i)**(-2))+x(2)*(V(i)**(-4/3))+x(3)*(V(i)**(-2/3))+x(4) end do open(unit=38,file="commment",status="replace" do i=1,n write(38,*)V(i),E(i),Ecal(i) end do write(38,*)V0,V0_2 close(38) stop end |
回复此楼» 本帖已获得的红花（最新10朵）
641960122» 猜你喜欢
今天务委会开完了，明天出结果吗
已经有17人回复
-
明天应该可查了！？
已经有6人回复
-
有没有大神帮我看看基金代码
已经有10人回复
-
放榜前的不淡定
已经有19人回复
-
某些机构，以效率低为荣，以效率低作为存在感
已经有10人回复
-
2026国自然函评费到账
已经有23人回复
-
在坚冰还盖着北海的时候，我看到了怒放的梅花。
已经有8人回复
-
今日不放榜？网传国自然预计 8 月 27 日可查结果
已经有22人回复
-
什么时候开奖？
已经有13人回复
-
如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于过柱子气泡的问题，请教前辈们
已经有17人回复
- 红外结果问题
已经有6人回复
- 新手无脑问题，希望有人解答~谢谢！（句法错误）
已经有6人回复
- 含有苄溴的硝基苯的还原问题
已经有10人回复
- APL 版面费的问题
已经有13人回复
- Chemical Reviews最新综述：水热/溶剂热直接溶胶凝胶法合成金属氧化物纳米结构
已经有792人回复
- 虚心请教！退火态是什么状态啊？
已经有8人回复
- 提供一个 Birch-Murnaghan 拟合的算法（更新1）
已经有95人回复
- 求绳子的曲线方程
已经有6人回复
- 关于声子谱对应的体积问题
已经有3人回复
- 求助椭偏拟合
已经有3人回复
- Matlab求解偏微分方程组
已经有22人回复
- 利用什么方法可以确定掺杂元素究竟进入晶格里面没有？
已经有18人回复
- 奇怪的投稿状态 International Journal of Computational Intelligence Systems
已经有14人回复
- 【求助】E-V曲线晶体状态方程拟合问题
已经有8人回复
- 【求助】如何从Murnaghan方程的结果中得到优化的晶格参数？
已经有11人回复
Gina88
木虫 (正式写手)
- 应助: 9 (幼儿园)
- 金币: 3328.8
- 散金: 124
- 红花: 8
- 帖子: 395
- 在线: 113小时
- 虫号: 1047684
- 注册: 2010-06-26
- 专业: 理论和计算化学
4楼2011-05-21 22:54:24
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
2楼2011-05-21 14:20:19
3楼2011-05-21 21:03:03
Gina88
木虫 (正式写手)
- 应助: 9 (幼儿园)
- 金币: 3328.8
- 散金: 124
- 红花: 8
- 帖子: 395
- 在线: 113小时
- 虫号: 1047684
- 注册: 2010-06-26
- 专业: 理论和计算化学
5楼2011-05-21 22:58:45
28