版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

当前主题已经存档。

zxzj05

荣誉版主 (著名写手)

应助: 0 (幼儿园)
贵宾: 1.242
金币: 5912.5
红花: 2
帖子: 1528
在线: 144.5小时
虫号: 406147
注册: 2007-06-18
性别: GG
专业: 金属功能材料
管辖: 第一性原理

[交流] 【分享】Code for creating EAM potential

EAM势函数是金属MD模拟非常有用的势函数，
其创建方法多种多样，
现在贴出一个C语言程序，和大家分享：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
/*
---------------------------------------------------------------
* Create Cu EAM potential file from Phys. Rev. B 63, 224106
* SMD, XMD and IMD format
*
* Compilation: g++ -DXXX -O mishincu.cpp -o mishincu
* Here XXX=SMD/XMD/IMD
*
* Usage: ./mishincu
* Then enter an inter number for the table size when prompted
* Copyright Sunny Bird
* sunzhmd@hotmail.com
---------------------------------------------------------------
*/
#include
#include
#include
#include
/*
==============================================
Fitting parameter for Cu potential functions
==============================================
*/
#define rc  5.50679
#define h  0.50037
#define E1  201.458
#define E2  6.59288e-3
#define r01  0.83591
#define r02  4.46867
#define r03  -2.19885
#define r04  -2.61984e2
#define afa1 2.97758
#define afa2 1.54927
#define delta 0.86225e-2
#define rs1  2.24000
#define rs2  1.80000
#define rs3  1.20000
#define S1  4.00000
#define S2  40.00000
#define S3  1.15000e3
#define a  3.80362
#define beta1 0.17394
#define beta2 5.35661e2
#define F0  -2.28235
#define F2  1.35535
#define q1  -1.27775
#define q2  -0.86074
#define q3  1.78804
#define q4  2.97571
#define Q1  0.40000
#define Q2  0.30000

/*
==================================
constants for differential function
==================================
*/
#define CON 1.4
#define CON2 (CON*CON)
#define BIG 1.0e30
#define NTAB 200
#define SAFE 2.0
double matrix[NTAB][NTAB];
/*
=========================
Mcros
=========================
*/
#define SEP ", \t\n"
#define FMAX(x,y) ((x>y)?(x)  y))
/*
===========================
Local subroutine prototype
===========================
*/
static double dfridr(double (*func)(double), double, double, double *);
static double pairfun(double);
static double denfun(double);
static double embfun(double);
int main( )
{
FILE *fp;
int n, nTable;
#if defined(SMD) || defined(XMD)
double rMin;
double deltaR;
#elif defined(IMD)
double r2Min, r2c;
double deltaR2;
double rr;
#endif
double rhoMax, deltaRho;
double r, rho;
double fpair, fden, femb;
#ifdef SMD
double dpair, dden, demb;
double errpair, errden, erremb;
#endif
/*
=========================================
Set parameter for tables
=========================================
*/
fprintf(stdout, " >>> Please input the size of the potential table to be created:");
scanf("%d", &nTable);
fprintf(stdout, " -  Thank you!\n\n");

#if defined(SMD) || defined(XMD)
rMin=1.0;
deltaR=(rc-rMin)/((double)(nTable-1));
#elif defined(IMD)
r2Min=1.0;
r2c=rc*rc;
deltaR2=(r2c-r2Min)/((double)(nTable-1));
#endif
rhoMax=2.0;
deltaRho=rhoMax/((double)(nTable-1));

#if defined(SMD) || defined(XMD)
if( (fp=fopen("cu.pot", "w"))==NULL)
{
  fprintf(stderr, " ***ERROR: Cannot create potential file cu.pot!***\n");
  fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
  exit(1);
}
#elif defined(IMD)
if( (fp=fopen("cu.pair", "w"))==NULL)
{
  fprintf(stderr, " ***ERROR: Cannot create pair potential file in IMD format!***\n");
  fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
  exit(1);
}
#else
fprintf(stderr, " ***ERROR: Please specify potential file format! ***\n");
fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
exit(1);
#endif
/*
===============================
Potential file header
===============================
*/
#if defined(SMD)
fprintf(fp, " -------------------------------------------\n");
fprintf(fp, " Cu EAM potential, Phys. Rev. B 63, 224106\n");
fprintf(fp, " -------------------------------------------\n");
fprintf(fp, " %5d %16.8e %16.8e %16.8e\n", nTable, rMin, rc, rhoMax);
#elif defined(XMD)
fprintf(fp, "#  -------------------------------------------\n");
fprintf(fp, "#  Cu EAM potential, Phys. Rev. B 63, 224106\n");
fprintf(fp, "#  -------------------------------------------\n");
fprintf(fp, "potential set eam 1\n");
fprintf(fp, "eunit ev\n");
fprintf(fp, "POTENTIAL PAIR 1 1 %d %16.8e %16.8e\n", nTable, rMin, rc);
#elif defined(IMD)
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "##  Cu EAM potential, Phys. Rev. B 63, 224106\n");
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "#F 2 1\n");
fprintf(fp, "#E\n");
fprintf(fp, "%16.8e %16.8e %16.8e\n", r2Min, r2c, deltaR2);
#endif
/*
=============================
Create pair potential table
=============================
*/
#if defined(SMD) || defined(XMD)
r=rMin;
#elif defined(IMD)
rr=r2Min;
#endif
for(n=0; n {
#if defined(SMD) || defined(XMD)
  fpair=pairfun(r);
#elif defined(IMD)
  r=sqrt(rr);
  fpair=pairfun(r);
#endif
#if defined(SMD)
  dpair=dfridr(pairfun, r, deltaR*5, &errpair);
  fden=denfun(r);
  dden=dfridr(denfun, r, deltaR*5, &errden);
  if(fabs(errpair)/fabs(dpair)<1e-0 && fabs(errden)/fabs(dden)<1e-0)
  {
fprintf(fp, "%16.8e %16.8e %16.8e %16.8e\n",
fpair, fden, dpair/r, dden/r); /* !!! This line defines the format of the potential file !!! */
  }
  else
  {
fprintf(stderr, " ***ERROR: Too big relative error for pairwise functions!***\n");
fprintf(stderr, " %10.7f %16.8e %16.8e %16.8e %16.8e\n", r,
fabs(errpair)/fabs(dpair), fabs(errden)/fabs(dden), dpair, dden);
fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
exit(1);
  }
#elif defined(XMD) || defined(IMD)
  fprintf(fp, "%16.8e\n",fpair);
#endif
#if defined(SMD) || defined(XMD)
  r+=deltaR;
#elif defined(IMD)
  rr+=deltaR2;
#endif
}
/*
================================
Create electron density table
================================
*/
#if defined(XMD)
fprintf(fp, "POTENTIAL DENS 1 %d %16.8e %16.8e\n", nTable, rMin, rc);
#elif defined(IMD)
fclose(fp);
if( (fp=fopen("cu.den", "w"))==NULL)
{
  fprintf(stderr, " ***ERROR: Cannot create electron density file in IMD format!***\n");
  fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
  exit(1);
}
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "##  Cu EAM potential, Phys. Rev. B 63, 224106\n");
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "#F 2 1\n");
fprintf(fp, "#E\n");
fprintf(fp, "%16.8e %16.8e %16.8e\n",r2Min, r2c,deltaR2);
#endif
#if defined(XMD)
r=rMin;
#elif defined(IMD)
rr=r2Min;
#endif

for(n=0; n {
#if defined(XMD)
  fden=denfun(r);
#elif defined(IMD)
  r=sqrt(rr);
  fden=denfun(r);
#endif
  fprintf(fp,"%16.8e\n",fden);
#if defined(XMD)
  r+=deltaR;
#elif defined(IMD)
  rr+=deltaR2;
#endif
}

/*
===============================
Create embedding energy table
===============================
*/
#if defined(XMD)
fprintf(fp, "POTENTIAL EMBED 1 %d %16.8e %16.8e\n", nTable, 0.0, rhoMax);
#elif defined(IMD)
fclose(fp);
if( (fp=fopen("cu.emb", "w"))==NULL)
{
  fprintf(stderr, " ***ERROR: Cannot create embedding energy file in IMD format!***\n");
  fprintf(stderr, " ---Location is line %d in file %s\n", __LINE__, __FILE__);
  exit(1);
}
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "##  Cu EAM potential, Phys. Rev. B 63, 224106\n");
fprintf(fp, "##  -------------------------------------------\n");
fprintf(fp, "#F 2 1\n");
fprintf(fp, "#E\n");
fprintf(fp, "%16.8e %16.8e %16.8e\n",0.0, rhoMax, deltaRho);
#endif
rho=0.0;
for(n=0; n {
  femb=embfun(rho);
#if defined(SMD)
  demb=dfridr(embfun, rho, deltaRho*10, &erremb);
  if(fabs(erremb)/fabs(demb)<1e-0)
  {
fprintf(fp, "%16.8e %16.8e\n", femb, demb);/* !!! This line defines the format of the potential file !!! */
  }
  else
  {
printf(" ***ERROR: Too big relative error for eam function!***\n");
printf(" %10.7f %16.8e %16.8e\n", rho, fabs(erremb)/fabs(demb), demb);
exit(1);
  }
#elif defined(XMD) || defined (IMD)
  fprintf(fp, "%16.8e\n",femb);
#endif
  rho+=deltaRho;
}
fclose(fp);
return 0;
}
/*
===============================
Subroutine to calculate
numeric direvatie
===============================
*/
double dfridr(double (*func)(double), double x, double dx, double *err)
{
int i,j;
double errt,fac,ans;
if (dx == 0.0)
{
  printf(" ***ERROR: dx must be non-zero!***\n");
  exit(1);
}
matrix[1][1]=((*func)(x+dx)-(*func)(x-dx))/(2.0*dx);
*err=BIG;
for (i=2;i<=NTAB;i++) {
  dx /= CON;
  matrix[1]=((*func)(x+dx)-(*func)(x-dx))/(2.0*dx);
  fac=CON2;
  for (j=2;j<=i;j++) {
matrix[j]=(matrix[j-1]*fac-matrix[j-1][i-1])/(fac-1.0);
fac=CON2*fac;
errt=FMAX(fabs(matrix[j]-matrix[j-1]),fabs(matrix[j]-matrix[j-1][i-1]));
if (errt <= *err) {
*err=errt;
ans=matrix[j];
}
  }
  if (fabs(matrix-matrix[i-1][i-1]) >= SAFE*(*err))
  {
break;
  }
}
return ans;
}
/*
=================================
Cu Pair Potential function
=================================
*/
static double pairfun(double r)
{
double morse1, morse2;
double psi;
double tempair;
double x;
morse1=exp(-2.0*afa1*(r-r01))-2.0*exp(-afa1*(r-r01));
morse2=exp(-2.0*afa2*(r-r02))-2.0*exp(-afa2*(r-r02));
x=(r-rc)/h;
if(x<0)
{
  psi=pow(x,4.0)/(1.0+pow(x,4.0));
}
else
{
  psi=0;
}
if(r {
  tempair=S1*pow(rs1-r,4.0)+S2*pow(rs2-r,4.0)+S3*pow(rs3-r,4.0);
}
else if(r {
  tempair=S1*pow(rs1-r,4.0)+S2*pow(rs2-r,4.0);
}
else if(r {
  tempair=S1*pow(rs1-r,4.0);
}
else
{
  tempair=0;
}
tempair=-tempair+(E1*morse1+E2*morse2+delta)*psi;
return tempair;
}
/*
===============================
Cu Electronic Dencity Function
===============================
*/
static double denfun(double r)
{
double tmpden;
double psi;
double x;
tmpden=a*exp(-beta1*(r-r03)*(r-r03))+exp(-beta2*(r-r04));
x=(r-rc)/h;
if(x<0)
{
  psi=pow(x,4.0)/(1.0+pow(x,4.0));
}
else
{
  psi=0;
}
tmpden*=psi;
return tmpden;
}
/*
====================================
Cu Embedding Function
====================================
*/
static double embfun(double den)
{
double tmpemb;;
tmpemb=F0+0.5*F2*(den-1.0)*(den-1.0);
if(den<=1.0)
{
  tmpemb=tmpemb
+q1*pow(den-1.0,3.0)
+q2*pow(den-1.0,4.0)
+q3*pow(den-1.0,5.0)
+q4*pow(den-1.0,6.0);
}
else
{
  tmpemb=tmpemb+q1*pow(den-1.0,3.0)+Q1*pow(den-1.0,4.0);
  tmpemb=tmpemb/(1.0+Q2*pow(den-1.0,3.0));
}
return tmpemb;
}

[ Last edited by lei0736 on 2009-11-25 at 13:49 ]

回复此楼

» 猜你喜欢

心脉受损已经有5人回复
博士读完未来一定会好吗已经有15人回复
Springer期刊投稿求助已经有4人回复
读博已经有3人回复
小论文投稿已经有3人回复
Bioresource Technology期刊，第一次返修的时候被退回好几次了已经有9人回复
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗已经有8人回复
申请2026年博士已经有6人回复

储氢家族欢迎储氢研究者！

1楼2009-06-04 15:28:37

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

043114076

木虫 (正式写手)

应助: 49 (小学生)
金币: 3364.9
散金: 52
红花: 23
帖子: 788
在线: 369.5小时
虫号: 570592
注册: 2008-06-07
性别: GG
专业: 理论和计算化学

好好研究一下，3q

回复此楼

2楼2009-06-14 00:23:00

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wuchenwf

荣誉版主 (职业作家)

应助: 0 (幼儿园)
贵宾: 3.433
金币: 19419.2
散金: 10
红花: 4
帖子: 3560
在线: 1035.7小时
虫号: 398569
注册: 2007-06-10
性别: GG
专业: 凝聚态物性I:结构、力学和
管辖: 第一性原理

★
小木虫(金币+0.5):给个红包，谢谢回帖交流

这么好的资源咋没人顶呢呵呵

赞一下(1人)

回复此楼

3楼2009-06-15 23:46:27

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wjhu

4楼2009-08-27 16:58:30

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

cxz314

木虫 (正式写手)

应助: 1 (幼儿园)
贵宾: 0.2
金币: 3815.1
散金: 993
红花: 1
帖子: 997
在线: 189小时
虫号: 475745
注册: 2007-12-07
性别: GG
专业: 凝聚态物性I:结构、力学和

★
小木虫(金币+0.5):给个红包，谢谢回帖交流

引用回帖:

Originally posted by zxzj05 at 2009-6-4 15:28:
EAM势函数是金属MD模拟非常有用的势函数，
其创建方法多种多样，
现在贴出一个C语言程序，和大家分享：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
/*
--------------------- ...

请问楼主，合金的EAM势怎么生成？

赞一下(1人)

回复此楼

5楼2009-08-28 11:03:56

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖