221 lines
7.1 KiB
C++
221 lines
7.1 KiB
C++
/* ----------------------------------------------------------------------
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
http://lammps.sandia.gov, Sandia National Laboratories
|
|
Steve Plimpton, sjplimp@sandia.gov
|
|
|
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
|
certain rights in this software. This software is distributed under
|
|
the GNU General Public License.
|
|
|
|
See the README file in the top-level LAMMPS directory.
|
|
------------------------------------------------------------------------- */
|
|
|
|
/* ----------------------------------------------------------------------
|
|
Contributing author: Aidan Thompson (SNL) - original Tersoff implementation
|
|
David Farrell (NWU) - ZBL addition
|
|
------------------------------------------------------------------------- */
|
|
|
|
#include "math.h"
|
|
#include "stdio.h"
|
|
#include "stdlib.h"
|
|
#include "string.h"
|
|
#include "pair_tersoff_zbl_cuda.h"
|
|
#include "atom.h"
|
|
#include "update.h"
|
|
#include "neighbor.h"
|
|
#include "neigh_list.h"
|
|
#include "neigh_request.h"
|
|
#include "force.h"
|
|
#include "comm.h"
|
|
#include "memory.h"
|
|
#include "error.h"
|
|
#include "math_const.h"
|
|
|
|
using namespace LAMMPS_NS;
|
|
using namespace MathConst;
|
|
|
|
#define MAXLINE 1024
|
|
#define DELTA 4
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
PairTersoffZBLCuda::PairTersoffZBLCuda(LAMMPS *lmp) : PairTersoffCuda(lmp)
|
|
{
|
|
// hard-wired constants in metal or real units
|
|
// a0 = Bohr radius
|
|
// epsilon0 = permittivity of vacuum = q / energy-distance units
|
|
// e = unit charge
|
|
// 1 Kcal/mole = 0.043365121 eV
|
|
|
|
if (strcmp(update->unit_style,"metal") == 0) {
|
|
global_a_0 = 0.529;
|
|
global_epsilon_0 = 0.00552635;
|
|
global_e = 1.0;
|
|
} else if (strcmp(update->unit_style,"real") == 0) {
|
|
global_a_0 = 0.529;
|
|
global_epsilon_0 = 0.00552635 * 0.043365121;
|
|
global_e = 1.0;
|
|
} else error->all(FLERR,"Pair tersoff/zbl requires metal or real units");
|
|
iszbl = true;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void PairTersoffZBLCuda::read_file(char *file)
|
|
{
|
|
int params_per_line = 21;
|
|
char **words = new char*[params_per_line+1];
|
|
|
|
delete [] params;
|
|
params = NULL;
|
|
nparams = 0;
|
|
|
|
// open file on proc 0
|
|
|
|
FILE *fp;
|
|
if (comm->me == 0) {
|
|
fp = fopen(file,"r");
|
|
if (fp == NULL) {
|
|
char str[128];
|
|
sprintf(str,"Cannot open Tersoff potential file %s",file);
|
|
error->one(FLERR,str);
|
|
}
|
|
}
|
|
|
|
// read each line out of file, skipping blank lines or leading '#'
|
|
// store line of params if all 3 element tags are in element list
|
|
|
|
int n,nwords,ielement,jelement,kelement;
|
|
char line[MAXLINE],*ptr;
|
|
int eof = 0;
|
|
|
|
while (1) {
|
|
if (comm->me == 0) {
|
|
ptr = fgets(line,MAXLINE,fp);
|
|
if (ptr == NULL) {
|
|
eof = 1;
|
|
fclose(fp);
|
|
} else n = strlen(line) + 1;
|
|
}
|
|
MPI_Bcast(&eof,1,MPI_INT,0,world);
|
|
if (eof) break;
|
|
MPI_Bcast(&n,1,MPI_INT,0,world);
|
|
MPI_Bcast(line,n,MPI_CHAR,0,world);
|
|
|
|
// strip comment, skip line if blank
|
|
|
|
if (ptr = strchr(line,'#')) *ptr = '\0';
|
|
nwords = atom->count_words(line);
|
|
if (nwords == 0) continue;
|
|
|
|
// concatenate additional lines until have params_per_line words
|
|
|
|
while (nwords < params_per_line) {
|
|
n = strlen(line);
|
|
if (comm->me == 0) {
|
|
ptr = fgets(&line[n],MAXLINE-n,fp);
|
|
if (ptr == NULL) {
|
|
eof = 1;
|
|
fclose(fp);
|
|
} else n = strlen(line) + 1;
|
|
}
|
|
MPI_Bcast(&eof,1,MPI_INT,0,world);
|
|
if (eof) break;
|
|
MPI_Bcast(&n,1,MPI_INT,0,world);
|
|
MPI_Bcast(line,n,MPI_CHAR,0,world);
|
|
if (ptr = strchr(line,'#')) *ptr = '\0';
|
|
nwords = atom->count_words(line);
|
|
}
|
|
|
|
if (nwords != params_per_line)
|
|
error->all(FLERR,"Incorrect format in Tersoff potential file");
|
|
|
|
// words = ptrs to all words in line
|
|
|
|
nwords = 0;
|
|
words[nwords++] = strtok(line," \t\n\r\f");
|
|
while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
|
|
|
|
// ielement,jelement,kelement = 1st args
|
|
// if all 3 args are in element list, then parse this line
|
|
// else skip to next line
|
|
|
|
for (ielement = 0; ielement < nelements; ielement++)
|
|
if (strcmp(words[0],elements[ielement]) == 0) break;
|
|
if (ielement == nelements) continue;
|
|
for (jelement = 0; jelement < nelements; jelement++)
|
|
if (strcmp(words[1],elements[jelement]) == 0) break;
|
|
if (jelement == nelements) continue;
|
|
for (kelement = 0; kelement < nelements; kelement++)
|
|
if (strcmp(words[2],elements[kelement]) == 0) break;
|
|
if (kelement == nelements) continue;
|
|
|
|
// load up parameter settings and error check their values
|
|
|
|
if (nparams == maxparam) {
|
|
maxparam += DELTA;
|
|
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
|
|
"pair:params");
|
|
}
|
|
|
|
params[nparams].ielement = ielement;
|
|
params[nparams].jelement = jelement;
|
|
params[nparams].kelement = kelement;
|
|
params[nparams].powerm = atof(words[3]);
|
|
params[nparams].gamma = atof(words[4]);
|
|
params[nparams].lam3 = atof(words[5]);
|
|
params[nparams].c = atof(words[6]);
|
|
params[nparams].d = atof(words[7]);
|
|
params[nparams].h = atof(words[8]);
|
|
params[nparams].powern = atof(words[9]);
|
|
params[nparams].beta = atof(words[10]);
|
|
params[nparams].lam2 = atof(words[11]);
|
|
params[nparams].bigb = atof(words[12]);
|
|
params[nparams].bigr = atof(words[13]);
|
|
params[nparams].bigd = atof(words[14]);
|
|
params[nparams].lam1 = atof(words[15]);
|
|
params[nparams].biga = atof(words[16]);
|
|
params[nparams].Z_i = atof(words[17]);
|
|
params[nparams].Z_j = atof(words[18]);
|
|
params[nparams].ZBLcut = atof(words[19]);
|
|
params[nparams].ZBLexpscale = atof(words[20]);
|
|
|
|
// currently only allow m exponent of 1 or 3
|
|
|
|
params[nparams].powermint = int(params[nparams].powerm);
|
|
|
|
if (
|
|
params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 ||
|
|
params[nparams].d < 0.0 || params[nparams].powern < 0.0 ||
|
|
params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
|
|
params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
|
|
params[nparams].bigd < 0.0 ||
|
|
params[nparams].bigd > params[nparams].bigr ||
|
|
params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
|
|
params[nparams].powerm - params[nparams].powermint != 0.0 ||
|
|
(params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
|
|
params[nparams].gamma < 0.0 ||
|
|
params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 ||
|
|
params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0)
|
|
error->all(FLERR,"Illegal Tersoff parameter");
|
|
|
|
nparams++;
|
|
}
|
|
|
|
delete [] words;
|
|
}
|
|
|
|
void PairTersoffZBLCuda::coeff(int narg, char **arg)
|
|
{
|
|
PairTersoffCuda::coeff(narg, arg);
|
|
for(int i=0;i<maxparam;i++)
|
|
{
|
|
params_f[i].a_ij = (0.8854*global_a_0) /
|
|
(pow(params[i].Z_i,0.23) + pow(params[i].Z_j,0.23));
|
|
params_f[i].premult = (params[i].Z_i * params[i].Z_j * pow(global_e,2.0))/(4.0*MY_PI*global_epsilon_0);
|
|
params_f[i].ZBLcut = params[i].ZBLcut;
|
|
params_f[i].ZBLexpscale = params[i].ZBLexpscale;
|
|
}
|
|
}
|