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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » Brich-Murnaghan状态方程去拟合平衡时的晶格体积出问题！
81/1返回列表
查看: 3177  |  回复: 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一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急 +3 Qj98d49dHFb9 2026-08-25 4/200 2026-08-26 03:08 by cNXvBfCpiZOM
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +17 木水思豆 2026-08-25 20/1000 2026-08-26 01:26 by chengyan1220
[基金申请] 今天务委会开完了，明天出结果吗 +18 angus9576 2026-08-25 22/1100 2026-08-26 00:41 by merchancy
[基金申请] 在坚冰还盖着北海的时候，我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 有没有大神帮我看看基金代码 28+4 1234567wang 2026-08-24 10/500 2026-08-25 19:15 by lfy8008
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 19/950 2026-08-25 17:57 by chick875
[基金申请] 2026国自然函评费到账 +20 羊腰板 2026-08-21 23/1150 2026-08-25 16:34 by zsna
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +5 process2012 2026-08-23 7/350 2026-08-25 09:42 by huixian257
[基金申请] 今天基金会出结果吗？20260819 +17 kkkl_v 2026-08-19 18/900 2026-08-25 09:41 by windflowerwy
[基金申请] 我面上完蛋了 +13 且听虎啸 2026-08-20 14/700 2026-08-25 09:10 by mrkang
[基金申请] 范进中举一文的中心思想 +7 炎黄贵胄 2026-08-22 8/400 2026-08-25 08:48 by ZJTJZ
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[基金申请] 今天放榜吗？ +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧？ +4 哈哈蛤？ 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
信息提示
关闭
请填处理意见
关闭