24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

牛排小黑娇

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 163.8
帖子: 49
在线: 34.5小时
虫号: 2728752
注册: 2013-10-16
性别: MM
专业: 金属非晶态、准晶和纳米晶

[求助] 势函数解析问题

我想做Mo-Ag的合金扩散模拟，但是并没有现成的势函数文件，我查阅文献查到了有人拟合过这个合金的势函数。

势函数形式是FS势格式和EAM势差不多的，我看lammps手册这个势函数可以通过一个DYNAMO的程序解析出来，不知道您是不是也做过这个势函数的解析文件，我想让您帮我看一下我的程序是不是有问题，我之前弄出来的模拟结果只能运行一行。希望尚老师您能抽时间帮我看下，万分感谢！

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <string>
#include <stdarg.h>
#include <string.h>
#include "vector"
using namespace std;
#define EOK 0
#define ERROR -1
double phi(double r,double cc,double c0,double c1,double c2,double c3,double c4)
{
double V_r;
if(r<=cc)
{
   V_r= (r-cc)*(r-cc)*(c0+c1*r+c2*r*r+c3*pow(r,3)+c4*pow(r,4));
}
else V_r=0.0;
return (V_r);
}
double rho(double r,double d,double beta)
{
double rho_r;
if(r<=d) rho_r= (r-d)*(r-d) + beta*beta*pow((r-d),4);
else rho_r=0.0;
return (rho_r);
}
double Frho(double rho,double AA)
{
double F_rho;
F_rho = -1.0*AA*pow(rho,0.5);
return (F_rho);
}
int main (void)
{
int Nr = 10000;
double rmax =9.8;                            这边取值的标准是？
double dr = rmax/(double)Nr;
int Nrho = 10000;
double rhomax = 45.5;
double drho = rhomax/(double)Nrho;
int i;
string  LAMMPSFilename="AgMo.eam.fs";
FILE *LAMMPSFile = fopen ((char*)LAMMPSFilename.data(), "w"

;
if (!LAMMPSFile) exit (ERROR);
// Header for setfl format
fprintf (LAMMPSFile, \
"#-> LAMMPS Potential File in DYNAMO 86 setfl Format <-#\n"\
"# Title :Embedded-Atom Method by Arthur F. Woter 1994\n"\
"# Implemented by Nize(2015) nize19890627@gmail.com\n"\
"2 Ag Mo\n"\
"%d %20.20f %d %20.20f %20.20f\n",Nrho,drho,Nr,dr,rmax);
fprintf (LAMMPSFile,"47 107.8682 4.0853 FCC\n"

;
// Embedding function and density function for Ag
   for(i=0;i<Nrho;i++)
   fprintf (LAMMPSFile,"%20.20e\n",Frho((double)i*drho,0.325514));
// Density function for Ag
  for (i = 0; i < Nr; i++)
fprintf (LAMMPSFile, "%20.20e\n",rho((double)i*dr,4.41,-1.293394));
fprintf (LAMMPSFile,"42 95.94 3.1472 BCC\n"

;
// Embedding function and density function for Mo
   for(i=0;i<Nrho;i++)
   fprintf (LAMMPSFile,"%20.20e\n",Frho((double)i*drho,1.848648));
// Density function for Mo
  for (i = 0; i < Nr; i++)
fprintf (LAMMPSFile, "%20.20e\n",rho((double)i*dr,4.14,0.0));
// Pair potential for Ag-Ag
  for (i = 0; i < Nr; i++)
      fprintf (LAMMPSFile, "%20.20e\n", (double)i*dr*phi((double)i*dr,4.76,10.681200,-12.045170,5.203072,-1.013304,0.0742308));
// Pair potential for Ag-Mo
  for (i = 0; i < Nr; i++)
      fprintf (LAMMPSFile, "%20.20e\n", (double)i*dr*phi((double)i*dr,4.50,44.406810,-45.490260,15.605110,-1.793704,0.0));
// Pair potential for Mo-Mo
  for (i = 0; i < Nr; i++)
      fprintf (LAMMPSFile, "%20.20e\n", (double)i*dr*phi((double)i*dr,3.2572,47.980660,-34.099240,5.832293,0.101749,0.0203934));
   fclose (LAMMPSFile);
// system("pause"

;
return (EOK);

知道的朋友帮我看看。急求谢谢啦