| 查看: 3186 | 回复: 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朵）
» 猜你喜欢
国自然评审意见
已经有10人回复
-
国自然面上复盘~欢迎讨论
已经有14人回复
-
系统查不到
已经有11人回复
-
有多少人是今天查系统知道结果的？
已经有18人回复
-
为什么资助数各大高校都创新高，自己申请怎么就这么难
已经有12人回复
-
有没有仍没收到信息的
已经有3人回复
-
导师吐槽：我怎么摊上了这么个极品研究生！
已经有8人回复
-
基金不中，共勉
已经有11人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有9人回复
-
2026年叶企孙基金
已经有6人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于过柱子气泡的问题，请教前辈们
已经有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人回复
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
- 专业: 理论和计算化学
4楼2011-05-21 22:54:24
Gina88
木虫 (正式写手)
- 应助: 9 (幼儿园)
- 金币: 3328.8
- 散金: 124
- 红花: 8
- 帖子: 395
- 在线: 113小时
- 虫号: 1047684
- 注册: 2010-06-26
- 专业: 理论和计算化学
5楼2011-05-21 22:58:45
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
6楼2011-05-21 23:03:56
Gina88
木虫 (正式写手)
- 应助: 9 (幼儿园)
- 金币: 3328.8
- 散金: 124
- 红花: 8
- 帖子: 395
- 在线: 113小时
- 虫号: 1047684
- 注册: 2010-06-26
- 专业: 理论和计算化学
|
program Brich_Murnaghan !announcement begin******************** implicit none real(8),allocatable::V( real(8),allocatable::E( real(8),allocatable::EV( real(8)::b(4),x(4) real(8)::a(4,4) real(8)::sum real(8)::E0,V0,B0,B0P real(8)::V0_2,B0_2,B0P_2 integer i,j,k,n,ios character(len=80)::nonsense real(8),allocatable::Ecal( !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,n b(i)=0;a(i,j)=0 end do do i=1,n a(2,1)=a(2,1)+V(i)**(-2/3) a(3,1)=a(3,1)+V(i)**(-4/3) a(4,1)=a(4,1)+V(i)**(-2) a(4,2)=a(4,2)+V(i)**(-8/3) a(4,3)=a(4,3)+V(i)**(-10/3) a(4,4)=a(4,4)+V(i)**(-4) b(1)=b(1)+E(i) b(2)=b(2)+E(i)*(V(i)**(-2/3)) b(3)=b(3)+E(i)*(V(i)**(-4/3)) b(4)=b(4)+E(i)*(V(i)**(-2)) end do a(1,1)=n a(1,2)=a(2,1) a(1,3)=a(3,1);a(2,2)=a(3,1) a(1,4)=a(4,1);a(2,3)=a(4,1);a(3,2)=a(4,1) a(2,4)=a(4,2);a(3,3)=a(4,2) a(3,4)=a(4,3) do k=1,3 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(1) sum=(42*x(4)*x(2)/(x(3)**2)-14)**2 sum=sum-(9*x(4)*x(2)/(x(3)**2)-3)*(196*x(4)*x(2)/(x(3)**2)-64) B0P=((42*x(4)*x(2)/(x(3)**2)-14)+sqrt(sum))/(9*x(4)*x(2)/(x(3)**2)-3) B0P_2=((42*x(4)*x(2)/(x(3)**2)-14)-sqrt(sum))/(9*x(4)*x(2)/(x(3)**2)-3) sum=-(x(4)/x(3))*((3*B0P-14)/(B0P-4)) if(sum>0)then V0=sum**1.5 B0=(16/9)*x(3)*(V0**(-7/3))/(14-3*B0P) end if sum=-(x(4)/x(3))*((3*B0P_2-14)/(B0P_2-4)) if (sum>0)then V0_2=sum**1.5 B0_2=(16/9)*x(3)*(V0_2**(-7/3))/(14-3*B0P_2) end if allocate(Ecal(n)) do i=1,n Ecal(i)=x(4)*(V(i)**(-2))+x(3)*(V(i)**(-4/3))+x(2)*(V(i)**(-2/3))+x(1) end do open(unit=38,file="comment",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 |
7楼2011-05-25 16:50:52
送红花一朵
 |
8楼2020-11-19 16:07:35
投票:
30