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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » LAMMPS/DL_POLY » 【分享】Code for creating EAM potential
71/1返回列表
查看: 1218  |  回复: 6
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

zxzj05

荣誉版主 (著名写手)

[交流] 【分享】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 ]
回复此楼

» 猜你喜欢
储氢家族欢迎储氢研究者！
1楼 2009-06-04 15:28:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

043114076

木虫 (正式写手)
好好研究一下，3q
回复此楼
2楼2009-06-14 00:23:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wuchenwf

荣誉版主 (职业作家)
★
小木虫(金币+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

木虫 (正式写手)
★
小木虫(金币+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的回帖

huangyc

木虫 (正式写手)
同问，期待解答
回复此楼
6楼2009-09-09 15:19:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bawdy

新虫 (小有名气)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
不错，顶，合金的咋整啊
赞一下(1人) 回复此楼
7楼2009-09-10 10:21:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zxzj05 的主题更新
71/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 一志愿985，本科211，0817化学工程与技术319求调剂 +5 Liwangman 2026-03-15 5/250 2026-03-16 17:10 by 我的船我的海
[考研] 0703化学336分求调剂 +3 zbzihdhd 2026-03-15 3/150 2026-03-16 16:44 by 我的船我的海
[考研] 0703一志愿211 285分求调剂 +5 ly3471z 2026-03-13 5/250 2026-03-16 16:16 by 哦哦123
[考研] 材料与化工求调剂 +3 为学666 2026-03-16 3/150 2026-03-16 15:09 by 加号+
[考研] 0703 物理化学调剂 +3 我可以上岸的对� 2026-03-13 5/250 2026-03-16 10:50 by 我可以上岸的对�
[考研] 22408总分284求调剂 +3 InAspic 2026-03-13 3/150 2026-03-15 11:10 by zhq0425
[考研] 289求调剂 +4 这么名字咋样 2026-03-14 6/300 2026-03-14 18:58 by userper
[考研] 297一志愿上交085600求调剂 +5 指尖八千里 2026-03-14 5/250 2026-03-14 17:26 by a不易
[考研] 297求调剂 +4 学海漂泊 2026-03-13 4/200 2026-03-14 11:51 by 热情沙漠
[考研] 331求调剂（0703有机化学 +5 ZY-05 2026-03-13 6/300 2026-03-14 10:51 by Jy?
[考研] 学硕285求调剂 +13 Wisjxn 2026-03-12 46/2300 2026-03-14 10:33 by JourneyLucky
[考研] 一志愿安徽大学材料工程专硕313分，求调剂的学校 +8 Yu先生 2026-03-10 10/500 2026-03-14 01:04 by JourneyLucky
[考研] 调剂 +3 13853210211 2026-03-10 3/150 2026-03-14 00:47 by JourneyLucky
[考研] 招收0805（材料）调剂 +3 18595523086 2026-03-13 3/150 2026-03-14 00:33 by 123%、
[考研] 311求调剂 +8 zchqwer 2026-03-10 8/400 2026-03-14 00:01 by JourneyLucky
[考研] 315求调剂 +9 小羊小羊_ 2026-03-11 10/500 2026-03-13 21:13 by SXNU李老师
[考研] 281求调剂 +9 Koxui 2026-03-12 11/550 2026-03-13 20:50 by Koxui
[考研] 293求调剂 +3 世界首富 2026-03-11 3/150 2026-03-13 16:27 by JourneyLucky
[考研] 290求调剂 +3 ADT 2026-03-13 3/150 2026-03-13 10:19 by peike
[考博] 26读博 +4 Rui135246 2026-03-12 10/500 2026-03-13 07:15 by gaobiao
信息提示
关闭
请填处理意见
关闭