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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » Brich-Murnaghan状态方程去拟合平衡时的晶格体积出问题！
81/1返回列表
查看: 3173  |  回复: 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返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 某些机构，以效率低为荣，以效率低作为存在感 +4 yuleib84 2026-08-25 5/250 2026-08-25 14:19 by xuqiangluo
[基金申请] 在坚冰还盖着北海的时候，我看到了怒放的梅花。 +3 ziyangfang 2026-08-25 3/150 2026-08-25 13:14 by huagongfeihu
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[基金申请] 2026国自然函评费到账 +19 羊腰板 2026-08-21 22/1100 2026-08-25 10:38 by popular289
[教师之家] 导师吐槽：我怎么摊上了这么个极品研究生！ +3 苏东坡二世 2026-08-23 3/150 2026-08-25 10:35 by shisan1313
[基金申请] 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
[基金申请] 明天应该可查了！？ +5 chengyan1220 2026-08-23 5/250 2026-08-24 23:28 by 我4大白菜
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[基金申请] filecode，4个jtjc了 +14 ziyangfang 2026-08-19 17/850 2026-08-24 18:37 by 哈哈蛤？
[基金申请] 让我中一个面上吧！ +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 17/850 2026-08-24 10:20 by echo8914667
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[教师之家] 跳槽后在研项目怎么办？ +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 今天放榜吗？ +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 时间戳今天，20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 应该是下周三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
信息提示
关闭
请填处理意见
关闭