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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » Brich-Murnaghan状态方程去拟合平衡时的晶格体积出问题！
81/1返回列表
查看: 2861  |  回复: 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返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 8rmuugja8q 2026-02-22 7/350 2026-02-23 09:44 by w4l55oybr1
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 khieu8v8m0 2026-02-22 8/400 2026-02-23 09:35 by w4l55oybr1
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 8/400 2026-02-23 09:29 by w4l55oybr1
[教师之家] 为什么中国大学工科教授们水了那么多所谓的顶会顶刊，但还是做不出宇树机器人？ +5 欢乐颂叶蓁 2026-02-21 8/400 2026-02-23 09:19 by 欢乐颂叶蓁
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 w89i99eaeh 2026-02-22 5/250 2026-02-23 08:04 by w4l55oybr1
[博后之家] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 6/300 2026-02-23 07:59 by w4l55oybr1
[博后之家] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +6 3dfhjxgsh7 2026-02-22 9/450 2026-02-23 07:49 by w4l55oybr1
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 4/200 2026-02-23 06:46 by jsjzfl
[硕博家园] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 8/400 2026-02-23 06:24 by w4l55oybr1
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +5 3dfhjxgsh7 2026-02-22 6/300 2026-02-23 02:04 by 5jlh3qtdvx
[教师之家] 版面费该交吗 +7 苹果在哪里 2026-02-22 8/400 2026-02-22 22:37 by otani
[基金申请] 面上可以超过30页吧？ +4 阿拉贡aragon 2026-02-22 4/200 2026-02-22 21:22 by 山西悬空寺空悬�
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 usprnugpzw 2026-02-21 6/300 2026-02-22 19:48 by w89i99eaeh
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 3dfhjxgsh7 2026-02-22 4/200 2026-02-22 16:52 by khieu8v8m0
[找工作] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 usprnugpzw 2026-02-22 3/150 2026-02-22 16:37 by khieu8v8m0
[公派出国] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 usprnugpzw 2026-02-21 4/200 2026-02-22 16:27 by khieu8v8m0
[基金申请] “人文社科而论，许多学术研究还没有达到民国时期的水平” +4 苏东坡二世 2026-02-18 5/250 2026-02-22 16:07 by liangep1573
[基金申请] 什么是人一生最重要的？ +4 瞬息宇宙 2026-02-21 4/200 2026-02-22 11:44 by huagongfeihu
[基金申请] 今年春晚有几个节目很不错，点赞！ +11 瞬息宇宙 2026-02-16 12/600 2026-02-21 21:14 by lq493392203
[基金申请] 体制内长辈说体制内绝大部分一辈子在底层，如同你们一样大部分普通教师忙且收入低 +9 瞬息宇宙 2026-02-20 12/600 2026-02-21 10:39 by 欢乐颂叶蓁
信息提示
关闭
请填处理意见
关闭