应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » Brich-Murnaghan状态方程去拟合平衡时的晶格体积出问题！
北京石油化工学院2026年研究生招生接收调剂公告
81/1返回列表
查看: 2911  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2011-05-21 00:02:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

【答案】应助回帖

★ ★
微尘、梦想(金币+2): 谢谢回复！ 2011-05-21 19:20:53
Gina88(金币+5): 那也谢谢！ 2011-05-21 22:51:40
这种问题其实没有办法回复，因为别人不可能看你给出的算法，再去核实你的程序……

我只看了一眼程序，觉得
CODE:
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

可能有点儿问题，因为你的两个分支的判断条件是一样的……
赞一下(1人) 回复此楼
2楼2011-05-21 14:20:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hookah

木虫 (正式写手)

【答案】应助回帖

★
Gina88(金币+2): 非常谢谢您 2011-05-21 22:52:17
余泽成(金币+1): 鼓励应助！ 2011-05-25 17:24:59
我于一个建议，作这个EOS的拟合，我建议你使用Origin或者MAtlab进行拟合，于EOS的方程了，不过是确定参数！
我作EOS就是用Origin进行拟合就行了！结果还可以！
赞一下(1人) 回复此楼
3楼2011-05-21 21:03:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Gina88

木虫 (正式写手)
引用回帖:
Originally posted by snoopyzhao at 2011-05-21 14:20:19:
这种问题其实没有办法回复，因为别人不可能看你给出的算法，再去核实你的程序……

我只看了一眼程序，觉得

[code]
if(sum>0)then
  V0=sum**1.5
  B0=(16/9)*x(2)*V0**(-7/3)/(14-3*B01)
  sum=-(x( ...

是不是算法都有一定的适用范围？
我用的是高斯消元法，
是不是EOS方程模拟晶体的E－V关系时，高斯消元法压根就不适用？
赞一下 回复此楼
4楼2011-05-21 22:54:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Gina88

木虫 (正式写手)
引用回帖:
Originally posted by snoopyzhao at 2011-05-21 14:20:19:
这种问题其实没有办法回复，因为别人不可能看你给出的算法，再去核实你的程序……

我只看了一眼程序，觉得

[code]
if(sum>0)then
  V0=sum**1.5
  B0=(16/9)*x(2)*V0**(-7/3)/(14-3*B01)
  sum=-(x( ...

我用的是附件中候柱峰的文档中EOS方程去拟合附件中的V-E关系，目的是得到平衡时的体积（能量最低时的体积）和体弹性模量。
算法是高斯消元法，最小二乘法拟合E－V数据。
谢谢大侠了
赞一下 回复此楼
5楼2011-05-21 22:58:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★
jjdg(金币+1): 辛苦了 2011-05-22 00:11:33
引用回帖:
Originally posted by Gina88 at 2011-05-21 22:58:45:
我用的是附件中候柱峰的文档中EOS方程去拟合附件中的V-E关系，目的是得到平衡时的体积（能量最低时的体积）和体弹性模量。
算法是高斯消元法，最小二乘法拟合E－V数据。
谢谢大侠了

俺不是你们那一行的，看不懂那个方程，实在抱歉……
赞一下(1人) 回复此楼
6楼2011-05-21 23:03:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Gina88

木虫 (正式写手)
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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

641960122
送红花一朵
8楼2020-11-19 16:07:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 Gina88 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料334求调剂 +18 Eecho# 2026-04-03 18/900 2026-04-06 00:56 by fmesaito
[考研] 332求调剂 +16 小小孟... 2026-04-05 17/850 2026-04-06 00:52 by fmesaito
[考研] 348求调剂 +6 wukira 2026-04-04 6/300 2026-04-05 18:11 by 猪会飞
[考研] 306分材料与化工求调剂 +7 黎吧啦啦你很有� 2026-04-03 7/350 2026-04-05 17:18 by Hdyxbekcb
[考研] 304求调剂（085602，过四级，一志愿985） +15 化工人999 2026-04-04 15/750 2026-04-05 16:25 by 我是电风扇r
[考研] 313求调剂 +5 海日海日 2026-04-04 7/350 2026-04-05 13:58 by imissbao
[考研] 271分求调剂学校 +12 zph158488！ 2026-04-02 13/650 2026-04-05 10:13 by lqwchd
[考研] 368求调剂 +5 今华习 2026-04-03 7/350 2026-04-04 18:47 by imissbao
[考研] 081200-11408-276学硕求调剂 +6 崔wj 2026-04-02 6/300 2026-04-03 10:19 by 蓝云思雨
[考研] 274求调剂 +10 薛定谔的虎。 2026-04-01 10/500 2026-04-03 10:13 by tianyyysss
[考研] 312 化工或制药调剂 +8 小小墨123 2026-04-02 9/450 2026-04-03 09:12 by zhouxiaoyu
[考研] 一志愿上海海洋大学083200食品学硕，求调剂，接受其他专业 +6 what张 2026-04-01 7/350 2026-04-02 16:48 by zzsw+
[考研] 311求调剂 +14 蓝月亮亮 2026-03-30 14/700 2026-04-02 12:18 by 1753564080
[考研] 0710生物学求调剂 +9 manman511 2026-04-01 9/450 2026-04-02 10:00 by zxl830724
[考研] 0710生物学，325求调剂 +3 mkkkkkl 2026-04-01 3/150 2026-04-02 09:48 by Jaylen.
[考研] 279求调剂 +6 学而思兮知 2026-04-01 6/300 2026-04-02 09:16 by vgtyfty
[考研] 考研调剂 +12 Amber00 2026-03-31 12/600 2026-04-02 09:04 by sanrepian
[考研] 一志愿西安交大材料学硕（英一数二）347，求调剂到高分子/材料相关专业 +7 zju51 2026-03-31 9/450 2026-04-01 19:35 by CFQZAFU
[考研] 349求调剂 +6 吃的不少 2026-04-01 6/300 2026-04-01 17:55 by JYD2011
[考研] 085601一志愿西北工业大学初试346 +4 085601初试346 2026-03-30 4/200 2026-03-31 07:47 by jp9609
信息提示
关闭
请填处理意见
关闭