jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6713 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-08-17 21:56:18 +00:00
parent 0ade235efa
commit 99c9091025
6 changed files with 2942 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1476
src/USER-MISC/fix_imd.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

131
src/USER-MISC/fix_imd.h Normal file
View File

@ -0,0 +1,131 @@
/* -*- 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
The FixIMD class contains code from VMD and NAMD which is copyrighted
by the Board of Trustees of the University of Illinois and is free to
use with LAMMPS according to point 2 of the UIUC license (10% clause):
" Licensee may, at its own expense, create and freely distribute
complimentary works that interoperate with the Software, directing others to
the TCBG server to license and obtain the Software itself. Licensee may, at
its own expense, modify the Software to make derivative works. Except as
explicitly provided below, this License shall apply to any derivative work
as it does to the original Software distributed by Illinois. Any derivative
work should be clearly marked and renamed to notify users that it is a
modified version and not the original Software distributed by Illinois.
Licensee agrees to reproduce the copyright notice and other proprietary
markings on any derivative work and to include in the documentation of such
work the acknowledgement:
"This software includes code developed by the Theoretical and Computational
Biophysics Group in the Beckman Institute for Advanced Science and
Technology at the University of Illinois at Urbana-Champaign."
Licensee may redistribute without restriction works with up to 1/2 of their
non-comment source code derived from at most 1/10 of the non-comment source
code developed by Illinois and contained in the Software, provided that the
above directions for notice and acknowledgement are observed. Any other
distribution of the Software or any derivative work requires a separate
license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to
negotiate an appropriate license for such distribution."
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(imd,FixIMD)
#else
#ifndef LMP_FIX_IMD_H
#define LMP_FIX_IMD_H
#include "fix.h"
#if defined(LAMMPS_ASYNC_IMD)
#include <pthread.h>
#endif
/* prototype for c wrapper that calls the real worker */
extern "C" void *fix_imd_ioworker(void *);
namespace LAMMPS_NS {
class FixIMD : public Fix {
public:
FixIMD(class LAMMPS *, int, char **);
~FixIMD();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
double memory_usage();
protected:
int imd_port;
void *localsock;
void *clientsock;
int num_coords; // total number of atoms controlled by this fix
int size_one; // bytes per atom in communication buffer.
int maxbuf; // size of atom communication buffer.
void *comm_buf; // communication buffer
void *idmap; // hash for mapping atom indices to consistent order.
int *rev_idmap; // list of the hash keys for reverse mapping.
int imd_forces; // number of forces communicated via IMD.
void *force_buf; // force data buffer
double imd_fscale; // scale factor for forces. in case VMD's units are off.
int imd_inactive; // true if IMD connection stopped.
int imd_terminate; // true if IMD requests termination of run.
int imd_trate; // IMD transmission rate.
int unwrap_flag; // true if coordinates need to be unwrapped before sending
int nowait_flag; // true if LAMMPS should not wait with the execution for VMD.
int connect_msg; // flag to indicate whether a "listen for connection message" is needed.
int me; // my MPI rank in this "world".
int nlevels_respa; // flag to determine respa levels.
int msglen;
char *msgdata;
#if defined(LAMMPS_ASYNC_IMD)
int buf_has_data; // flag to indicate to the i/o thread what to do.
pthread_mutex_t write_mutex; // mutex for sending coordinates to i/o thread
pthread_cond_t write_cond; // conditional variable for coordinate i/o
pthread_mutex_t read_mutex; // mutex for accessing data receieved by i/o thread
pthread_t iothread; // thread id for i/o thread.
pthread_attr_t iot_attr; // i/o thread attributes.
public:
void ioworker(void);
#endif
protected:
int reconnect();
};
}
#endif
#endif
// Local Variables:
// mode: c++
// compile-command: "make -j4 openmpi"
// c-basic-offset: 2
// fill-column: 76
// indent-tabs-mode: nil
// End:

403
src/USER-MISC/fix_smd.cpp Normal file
View File

@ -0,0 +1,403 @@
/* ----------------------------------------------------------------------
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: Axel Kohlmeyer (UPenn)
based on fix spring by: Paul Crozier (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_smd.h"
#include "atom.h"
#include "comm.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
#include "group.h"
using namespace LAMMPS_NS;
enum { SMD_NONE=0,
SMD_TETHER=1<<0, SMD_COUPLE=1<<1,
SMD_CVEL=1<<2, SMD_CFOR=1<<3,
SMD_AUTOX=1<<4, SMD_AUTOY=1<<5, SMD_AUTOZ=1<<6};
#define SMALL 0.001
/* ---------------------------------------------------------------------- */
FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
styleflag = SMD_NONE;
k_smd = f_smd = v_smd = -1.0;
xflag = yflag = zflag = 1;
xc = yc = zc = 0.0;
xn = yn = zn = 1.0;
pmf = r_old = r_now = r0 = 0.0;
restart_global = 1;
vector_flag = 1;
size_vector = 7;
global_freq = 1;
extvector = 1;
int argoffs=3;
if (strcmp(arg[argoffs],"cvel") == 0) {
if (narg < argoffs+3) error->all("Illegal fix smd command");
styleflag |= SMD_CVEL;
k_smd = atof(arg[argoffs+1]);
v_smd = atof(arg[argoffs+2]); // to be multiplied by update->dt when used.
argoffs += 3;
} else if (strcmp(arg[argoffs],"cfor") == 0) {
if (narg < argoffs+2) error->all("Illegal fix smd command");
styleflag |= SMD_CFOR;
f_smd = atof(arg[argoffs+1]);
argoffs += 2;
} else error->all("Illegal fix smd command");
if (strcmp(arg[argoffs],"tether") == 0) {
if (narg < argoffs+5) error->all("Illegal fix smd command");
styleflag |= SMD_TETHER;
if (strcmp(arg[argoffs+1],"NULL") == 0) xflag = 0;
else xc = atof(arg[argoffs+1]);
if (strcmp(arg[argoffs+2],"NULL") == 0) yflag = 0;
else yc = atof(arg[argoffs+2]);
if (strcmp(arg[argoffs+3],"NULL") == 0) zflag = 0;
else zc = atof(arg[argoffs+3]);
r0 = atof(arg[argoffs+4]);
if (r0 < 0) error->all("R0 < 0 for fix smd command");
argoffs += 5;
} else if (strcmp(arg[argoffs],"couple") == 0) {
if (narg < argoffs+6) error->all("Illegal fix smd command");
styleflag |= SMD_COUPLE;
igroup2 = group->find(arg[argoffs+1]);
if (igroup2 == -1)
error->all("Could not find fix smd couple group ID");
if (igroup2 == igroup)
error->all("Two groups cannot be the same in fix smd couple");
group2bit = group->bitmask[igroup2];
if (strcmp(arg[argoffs+2],"NULL") == 0) xflag = 0;
else if (strcmp(arg[argoffs+2],"auto") == 0) styleflag |= SMD_AUTOX;
else xc = atof(arg[argoffs+2]);
if (strcmp(arg[argoffs+3],"NULL") == 0) yflag = 0;
else if (strcmp(arg[argoffs+3],"auto") == 0) styleflag |= SMD_AUTOY;
else yc = atof(arg[argoffs+3]);
if (strcmp(arg[argoffs+4],"NULL") == 0) zflag = 0;
else if (strcmp(arg[argoffs+4],"auto") == 0) styleflag |= SMD_AUTOZ;
else zc = atof(arg[argoffs+4]);
r0 = atof(arg[argoffs+5]);
if (r0 < 0) error->all("R0 < 0 for fix smd command");
argoffs +=6;
} else error->all("Illegal fix smd command");
force_flag = 0;
ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
}
/* ---------------------------------------------------------------------- */
int FixSMD::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMD::init()
{
double xcm[3], xcm2[3];
masstotal = group->mass(igroup);
group->xcm(igroup,masstotal,xcm);
double dx,dy,dz;
if (styleflag & SMD_TETHER) {
dx = xc - xcm[0];
dy = yc - xcm[1];
dz = zc - xcm[2];
} else { /* SMD_COUPLE */
masstotal2 = group->mass(igroup2);
group->xcm(igroup2,masstotal2,xcm2);
if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0];
else dx = xc;
if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1];
else dy = yc;
if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2];
else dz = zc;
}
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;
r_old = sqrt(dx*dx + dy*dy + dz*dz);
if (r_old > SMALL) {
xn = dx/r_old;
yn = dy/r_old;
zn = dz/r_old;
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixSMD::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixSMD::post_force(int vflag)
{
if (styleflag & SMD_TETHER) smd_tether();
else smd_couple();
if (styleflag & SMD_CVEL) r_old += v_smd * update->dt;
}
/* ---------------------------------------------------------------------- */
void FixSMD::smd_tether()
{
double xcm[3];
group->xcm(igroup,masstotal,xcm);
// fx,fy,fz = components of k * (r-r0)
double dx,dy,dz,fx,fy,fz,r,dr;
dx = xcm[0] - xc;
dy = xcm[1] - yc;
dz = xcm[2] - zc;
r_now = sqrt(dx*dx + dy*dy + dz*dz);
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;
r = sqrt(dx*dx + dy*dy + dz*dz);
if (styleflag & SMD_CVEL) {
if(r > SMALL) {
double fsign;
fsign = (v_smd<0.0) ? -1.0 : 1.0;
dr = r - r0 - r_old;
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;
}
} else {
r_old = r;
fx = f_smd*dx/r;
fy = f_smd*dy/r;
fz = f_smd*dz/r;
}
// apply restoring force to atoms in group
// f = -k*(r-r0)*mass/masstotal
double **f = atom->f;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
int nlocal = atom->nlocal;
ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
force_flag = 0;
double massfrac;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massfrac = mass[type[i]]/masstotal;
f[i][0] -= fx*massfrac;
f[i][1] -= fy*massfrac;
f[i][2] -= fz*massfrac;
ftotal[0] -= fx*massfrac;
ftotal[1] -= fy*massfrac;
ftotal[2] -= fz*massfrac;
}
}
/* ---------------------------------------------------------------------- */
void FixSMD::smd_couple()
{
double xcm[3],xcm2[3];
group->xcm(igroup,masstotal,xcm);
group->xcm(igroup2,masstotal2,xcm2);
// renormalize direction of spring
double dx,dy,dz,r,dr;
if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0];
else dx = xn*r_old;
if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1];
else dy = yn*r_old;
if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2];
else dz = zn*r_old;
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;
r = sqrt(dx*dx + dy*dy + dz*dz);
if (r > SMALL) {
xn = dx/r; yn = dy/r; zn = dz/r;
}
double fx,fy,fz;
if (styleflag & SMD_CVEL) {
dx = xcm2[0] - xcm[0];
dy = xcm2[1] - xcm[1];
dz = xcm2[2] - xcm[2];
r_now = sqrt(dx*dx + dy*dy + dz*dz);
dx -= xn*r_old;
dy -= yn*r_old;
dz -= zn*r_old;
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;
r = sqrt(dx*dx + dy*dy + dz*dz);
dr = r - r0;
if (r > SMALL) {
double fsign;
fsign = (v_smd<0.0) ? -1.0 : 1.0;
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt;
}
} else {
dx = xcm2[0] - xcm[0];
dy = xcm2[1] - xcm[1];
dz = xcm2[2] - xcm[2];
r_now = sqrt(dx*dx + dy*dy + dz*dz);
r_old = r;
fx = f_smd*xn;
fy = f_smd*yn;
fz = f_smd*zn;
}
// apply restoring force to atoms in group
// f = -k*(r-r0)*mass/masstotal
double **f = atom->f;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
int nlocal = atom->nlocal;
ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
force_flag = 0;
double massfrac;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
massfrac = mass[type[i]]/masstotal;
f[i][0] += fx*massfrac;
f[i][1] += fy*massfrac;
f[i][2] += fz*massfrac;
ftotal[0] += fx*massfrac;
ftotal[1] += fy*massfrac;
ftotal[2] += fz*massfrac;
}
if (mask[i] & group2bit) {
massfrac = mass[type[i]]/masstotal2;
f[i][0] -= fx*massfrac;
f[i][1] -= fy*massfrac;
f[i][2] -= fz*massfrac;
}
}
}
/* ---------------------------------------------------------------------- */
void FixSMD::write_restart(FILE *fp)
{
#define RESTART_ITEMS 5
double buf[RESTART_ITEMS], fsign;
if (comm->me == 0) {
// make sure we project the force into the direction of the pulling.
fsign = (v_smd<0.0) ? -1.0 : 1.0;
buf[0] = r_old;
buf[1] = xn*fsign;
buf[2] = yn*fsign;
buf[3] = zn*fsign;
buf[4] = pmf;
int size = RESTART_ITEMS*sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(&buf[0],sizeof(double),RESTART_ITEMS,fp);
}
}
/* ---------------------------------------------------------------------- */
void FixSMD::restart(char *buf)
{
double *list = (double *)buf;
r_old = list[0];
xn=list[1];
yn=list[2];
zn=list[3];
pmf=list[4];
}
/* ---------------------------------------------------------------------- */
void FixSMD::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ----------------------------------------------------------------------
return components of total smd force on fix group
------------------------------------------------------------------------- */
double FixSMD::compute_vector(int n)
{
// only sum across procs one time
if (force_flag == 0) {
MPI_Allreduce(ftotal,ftotal_all,3,MPI_DOUBLE,MPI_SUM,world);
force_flag = 1;
if (styleflag & SMD_CVEL) {
ftotal_all[3]=ftotal_all[0]*xn+ftotal_all[1]*yn+ftotal_all[2]*zn;
ftotal_all[4]=r_old;
} else {
ftotal_all[3]=f_smd;
ftotal_all[4]=r_old;
}
ftotal_all[5]=r_now;
ftotal_all[6]=pmf;
}
return ftotal_all[n];
}

60
src/USER-MISC/fix_smd.h Normal file
View File

@ -0,0 +1,60 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd,FixSMD)
#else
#ifndef LMP_FIX_SMD_H
#define LMP_FIX_SMD_H
#include "fix.h"
namespace LAMMPS_NS {
class FixSMD : public Fix {
public:
FixSMD(class LAMMPS *, int, char **);
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
double compute_vector(int);
void write_restart(FILE *);
void restart(char *);
private:
double xc,yc,zc,xn,yn,zn,r0;
double k_smd,f_smd,v_smd;
int xflag,yflag,zflag;
int styleflag;
double r_old,r_now,pmf;
int igroup2,group2bit;
double masstotal,masstotal2;
int nlevels_respa;
double ftotal[3],ftotal_all[7];
int force_flag;
void smd_tether();
void smd_couple();
};
}
#endif
#endif

642
src/USER-MISC/pair_cdeam.cpp Normal file
View File

@ -0,0 +1,642 @@
/* ----------------------------------------------------------------------
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: Alexander Stukowski
Technical University of Darmstadt,
Germany Department of Materials Science
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_cdeam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
// This is for debugging purposes. The ASSERT() macro is used in the code to check
// if everything runs as expected. Change this to #if 0 if you don't need the checking.
#if 0
#define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop())
inline void my_noop() {}
inline void my_failure(Error* error, const char* file, int line) {
char str[1024];
sprintf(str,"Assertion failure: File %s, line %i", file, line);
error->one(str);
}
#else
#define ASSERT(cond)
#endif
#define MAXLINE 1024 // This sets the maximum line length in EAM input files.
PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)
{
single_enable = 0;
rhoB = NULL;
D_values = NULL;
hcoeff = NULL;
// Set communication buffer sizes needed by this pair style.
if(cdeamVersion == 1) {
comm_forward = 4;
comm_reverse = 3;
}
else if(cdeamVersion == 2) {
comm_forward = 3;
comm_reverse = 2;
}
else {
error->all("Invalid CD-EAM potential version.");
}
}
PairCDEAM::~PairCDEAM()
{
memory->destroy(rhoB);
memory->destroy(D_values);
if(hcoeff) delete[] hcoeff;
}
void PairCDEAM::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rhoip,rhojp,recip,phi;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
// Grow per-atom arrays if necessary
if(atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(fp);
memory->destroy(rhoB);
memory->destroy(D_values);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rhoB,nmax,"pair:rhoB");
memory->create(fp,nmax,"pair:fp");
memory->create(D_values,nmax,"pair:D_values");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Zero out per-atom arrays.
int m = nlocal + atom->nghost;
for(i = 0; i < m; i++) {
rho[i] = 0.0;
rhoB[i] = 0.0;
D_values[i] = 0.0;
}
// Stage I
// Compute rho and rhoB at each local atom site.
// Additionally calculate the D_i values here if we are using the one-site formulation.
// For the two-site formulation we have to calculate the D values in an extra loop (Stage II).
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for(jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
double localrho = RhoOfR(index, jtype, itype);
rho[i] += localrho;
if(jtype == speciesB) rhoB[i] += localrho;
if(newton_pair || j < nlocal) {
localrho = RhoOfR(index, itype, jtype);
rho[j] += localrho;
if(itype == speciesB) rhoB[j] += localrho;
}
if(cdeamVersion == 1 && itype != jtype) {
// Note: if the i-j interaction is not concentration dependent (because either
// i or j are not species A or B) then its contribution to D_i and D_j should
// be ignored.
// This if-clause is only required for a ternary.
if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) {
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
D_values[i] += Phi_AB;
if(newton_pair || j < nlocal)
D_values[j] += Phi_AB;
}
}
}
}
}
// Communicate and sum densities.
if(newton_pair) {
communicationStage = 1;
comm->reverse_comm_pair(this);
}
// fp = derivative of embedding energy at each atom
// phi = embedding energy at each atom
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
EAMTableIndex index = rhoToTableIndex(rho[i]);
fp[i] = FPrimeOfRho(index, type[i]);
if(eflag) {
phi = FofRho(index, type[i]);
if (eflag_global) eng_vdwl += phi;
if (eflag_atom) eatom[i] += phi;
}
}
// Communicate derivative of embedding function and densities
// and D_values (this for one-site formulation only).
communicationStage = 2;
comm->forward_comm_pair(this);
// The electron densities may not drop to zero because then the concentration would no longer be defined.
// But the concentration is not needed anyway if there is no interaction with another atom, which is the case
// if the electron density is exactly zero. That's why the following lines have been commented out.
//
//for(i = 0; i < nlocal + atom->nghost; i++) {
// if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
// error->one("CD-EAM potential routine: Detected atom with zero electron density.");
//}
// Stage II
// This is only required for the original two-site formulation of the CD-EAM potential.
if(cdeamVersion == 2) {
// Compute intermediate value D_i for each atom.
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// This code line is required for ternary alloys.
if(itype != speciesA && itype != speciesB) continue;
double x_i = rhoB[i] / rho[i]; // Concentration at atom i.
for(jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if(itype == jtype) continue;
// This code line is required for ternary alloys.
if(jtype != speciesA && jtype != speciesB) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// The concentration independent part of the cross pair potential.
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
// Average concentration of two sites
double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]);
// Calculate derivative of h(x_ij) polynomial function.
double h_prime = evalHprime(x_ij);
D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]);
if(newton_pair || j < nlocal)
D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
}
}
}
// Communicate and sum D values.
if(newton_pair) {
communicationStage = 3;
comm->reverse_comm_pair(this);
}
communicationStage = 4;
comm->forward_comm_pair(this);
}
// Stage III
// Compute force acting on each atom.
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// Concentration at site i
double x_i = -1.0; // The value -1 indicates: no concentration dependence for all interactions of atom i.
// It will be replaced by the concentration at site i if atom i is either A or B.
double D_i, h_prime_i;
// This if-clause is only required for ternary alloys.
if((itype == speciesA || itype == speciesB) && rho[i] != 0.0) {
// Compute local concentration at site i.
x_i = rhoB[i]/rho[i];
ASSERT(x_i >= 0 && x_i<=1.0);
if(cdeamVersion == 1) {
// Calculate derivative of h(x_i) polynomial function.
h_prime_i = evalHprime(x_i);
D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
}
else if(cdeamVersion == 2) {
D_i = D_values[i];
}
else ASSERT(false);
}
for(jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
// psip needs both fp[i] and fp[j] terms since r_ij appears in two
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
rhoip = RhoPrimeOfR(index, itype, jtype);
rhojp = RhoPrimeOfR(index, jtype, itype);
fpair = fp[i]*rhojp + fp[j]*rhoip;
recip = 1.0/r;
double x_j = -1; // The value -1 indicates: no concentration dependence for this i-j pair
// because atom j is not of species A nor B.
// This code line is required for ternary alloy.
if(jtype == speciesA || jtype == speciesB) {
ASSERT(rho[i] != 0.0);
ASSERT(rho[j] != 0.0);
// Compute local concentration at site j.
x_j = rhoB[j]/rho[j];
ASSERT(x_j >= 0 && x_j<=1.0);
double D_j;
if(cdeamVersion == 1) {
// Calculate derivative of h(x_j) polynomial function.
double h_prime_j = evalHprime(x_j);
D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
}
else if(cdeamVersion == 2) {
D_j = D_values[j];
}
else ASSERT(false);
double t2 = -rhoB[j];
if(itype == speciesB) t2 += rho[j];
fpair += D_j * rhoip * t2;
}
// This if-clause is only required for a ternary alloy.
// Actually we don't need it at all because D_i should be zero anyway if
// atom i has no concentration dependent interactions (because it is not species A or B).
if(x_i != -1.0) {
double t1 = -rhoB[i];
if(jtype == speciesB) t1 += rho[i];
fpair += D_i * rhojp * t1;
}
double phip;
double phi = PhiOfR(index, itype, jtype, recip, phip);
if(itype == jtype || x_i == -1.0 || x_j == -1.0) {
// Case of no concentration dependence.
fpair += phip;
}
else {
// We have a concentration dependence for the i-j interaction.
double h;
if(cdeamVersion == 1) {
// Calculate h(x_i) polynomial function.
double h_i = evalH(x_i);
// Calculate h(x_j) polynomial function.
double h_j = evalH(x_j);
h = 0.5 * (h_i + h_j);
}
else if(cdeamVersion == 2) {
// Average concentration.
double x_ij = 0.5 * (x_i + x_j);
// Calculate h(x_ij) polynomial function.
h = evalH(x_ij);
}
else ASSERT(false);
fpair += h * phip;
phi *= h;
}
// Divide by r_ij and negate to get forces from gradient.
fpair /= -r;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if(newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if(eflag) evdwl = phi;
if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if(vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::coeff(int narg, char **arg)
{
PairEAMAlloy::coeff(narg, arg);
// Make sure the EAM file is a CD-EAM binary alloy.
if(setfl->nelements < 2)
error->all("The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");
// Read in the coefficients of the h polynomial from the end of the EAM file.
read_h_coeff(arg[2]);
// Determine which atom type is the A species and which is the B species in the alloy.
// By default take the first element (index 0) in the EAM file as the A species
// and the second element (index 1) in the EAM file as the B species.
speciesA = -1;
speciesB = -1;
for(int i = 1; i <= atom->ntypes; i++) {
if(map[i] == 0) {
if(speciesA >= 0)
error->all("The first element from the EAM file may only be mapped to a single atom type.");
speciesA = i;
}
if(map[i] == 1) {
if(speciesB >= 0)
error->all("The second element from the EAM file may only be mapped to a single atom type.");
speciesB = i;
}
}
if(speciesA < 0)
error->all("The first element from the EAM file must be mapped to exactly one atom type.");
if(speciesB < 0)
error->all("The second element from the EAM file must be mapped to exactly one atom type.");
}
/* ----------------------------------------------------------------------
Reads in the h(x) polynomial coefficients
------------------------------------------------------------------------- */
void PairCDEAM::read_h_coeff(char *filename)
{
if(comm->me == 0) {
// Open potential file
FILE *fp;
char line[MAXLINE];
char nextline[MAXLINE];
fp = fopen(filename,"r");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open EAM potential file %s", filename);
error->one(str);
}
// h coefficients are stored at the end of the file.
// Skip to last line of file.
while(fgets(nextline, MAXLINE, fp) != NULL) {
strcpy(line, nextline);
}
char* ptr = strtok(line, " \t\n\r\f");
int degree = atoi(ptr);
nhcoeff = degree+1;
hcoeff = new double[nhcoeff];
int i = 0;
while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
hcoeff[i++] = atof(ptr);
}
if(i != nhcoeff || nhcoeff < 1)
error->one("Failed to read h(x) function coefficients from EAM file.");
// Close the potential file.
fclose(fp);
}
MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
if(comm->me != 0) hcoeff = new double[nhcoeff];
MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
}
/* ---------------------------------------------------------------------- */
int PairCDEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if(communicationStage == 2) {
if(cdeamVersion == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
buf[m++] = D_values[j];
}
return 4;
}
else if(cdeamVersion == 2) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
}
return 3;
}
else { ASSERT(false); return 0; }
}
else if(communicationStage == 4) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = D_values[j];
}
return 1;
}
else return 0;
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if(communicationStage == 2) {
if(cdeamVersion == 1) {
for(i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
D_values[i] = buf[m++];
}
}
else if(cdeamVersion == 2) {
for(i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
}
}
else ASSERT(false);
}
else if(communicationStage == 4) {
for(i = first; i < last; i++) {
D_values[i] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int PairCDEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if(communicationStage == 1) {
if(cdeamVersion == 1) {
for(i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
buf[m++] = D_values[i];
}
return 3;
}
else if(cdeamVersion == 2) {
for(i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
}
return 2;
}
else { ASSERT(false); return 0; }
}
else if(communicationStage == 3) {
for(i = first; i < last; i++) {
buf[m++] = D_values[i];
}
return 1;
}
else return 0;
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if(communicationStage == 1) {
if(cdeamVersion == 1) {
for(i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
D_values[j] += buf[m++];
}
}
else if(cdeamVersion == 2) {
for(i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
}
}
else ASSERT(false);
}
else if(communicationStage == 3) {
for(i = 0; i < n; i++) {
j = list[i];
D_values[j] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairCDEAM::memory_usage()
{
double bytes = 2 * nmax * sizeof(double);
return PairEAMAlloy::memory_usage() + bytes;
}

230
src/USER-MISC/pair_cdeam.h Normal file
View File

@ -0,0 +1,230 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(eam/cd,PairCDEAM_OneSite)
PairStyle(eam/cd/old,PairCDEAM_TwoSite)
#else
#ifndef LMP_PAIR_CDEAM_H
#define LMP_PAIR_CDEAM_H
#include "pair_eam_alloy.h"
namespace LAMMPS_NS {
class PairCDEAM : public PairEAMAlloy
{
public:
/// Constructor.
PairCDEAM(class LAMMPS*, int cdeamVersion);
/// Destructor.
virtual ~PairCDEAM();
/// Calculates the energies and forces for all atoms in the system.
void compute(int, int);
/// Parses the pair_coeff command parameters for this pair style.
void coeff(int, char **);
/// This is for MPI communication with neighbor nodes.
int pack_comm(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
/// Reports the memory usage of this pair style to LAMMPS.
double memory_usage();
/// Parses the coefficients of the h polynomial from the end of the EAM file.
void read_h_coeff(char* filename);
public:
// The public interface exposed by this potential class.
// Evaluates the h(x) polynomial for a given local concentration x.
inline double evalH(double x) const {
double v = 0.0;
for(int i = nhcoeff-1; i >= 1; i--) {
v = (v + hcoeff[i]) * x;
}
return v + hcoeff[0];
}
// Calculates the derivative of the h(x) polynomial.
inline double evalHprime(double x) const {
double v = 0.0;
for(int i = nhcoeff-1; i >= 2; i--) {
v = (v + (double)i * hcoeff[i]) * x;
}
return v + hcoeff[1];
}
// We have two versions of the CD-EAM potential. The user can choose which one he wants to use:
//
// Version 1 - One-site concentration: The h(x_i) function depends only on the concentration at the atomic site i.
// This is a new version with a slight modification to the formula. It happens to be computationally more efficient.
// It has been published in the MSMSE 2009 paper.
//
// Version 2 - Two-site concentration: The h(x_ij) function depends on the concentrations at two atomic sites i and j.
// This is the original version from the 2005 PRL paper.
int cdeamVersion;
// Per-atom arrays
// The partial density of B atoms at each atom site.
double *rhoB;
// The intermediate value D_i for each atom.
// The meaning of these values depend on the version of the CD-EAM potential used:
//
// For the one-site concentration version these are the v_i values defined in equation (21)
// of the MSMSE paper.
//
// For the old two-site concentration version these are the M_i values defined in equation (12)
// of the MSMSE paper.
double *D_values;
// The atom type index that is considered to be the A species in the alloy.
int speciesA;
// The atom type index that is considered to be the B species in the alloy.
int speciesB;
protected:
// Evaluation functions:
// This structure specifies an entry in one of the EAM spline tables
// and the corresponding floating point part.
typedef struct {
int m;
double p;
} EAMTableIndex;
// Converts a radius value to an index value to be used in a spline table lookup.
inline EAMTableIndex radiusToTableIndex(double r) const {
EAMTableIndex index;
index.p = r*rdr + 1.0;
index.m = static_cast<int>(index.p);
index.m = index.m <= (nr-1) ? index.m : (nr-1);
index.p -= index.m;
index.p = index.p <= 1.0 ? index.p : 1.0;
return index;
}
// Converts a density value to an index value to be used in a spline table lookup.
inline EAMTableIndex rhoToTableIndex(double rho) const {
EAMTableIndex index;
index.p = rho*rdrho + 1.0;
index.m = static_cast<int>(index.p);
index.m = index.m <= (nrho-1) ? index.m : (nrho-1);
index.p -= index.m;
index.p = index.p <= 1.0 ? index.p : 1.0;
return index;
}
// Computes the derivative of rho(r)
inline double RhoPrimeOfR(const EAMTableIndex& index, int itype, int jtype) const {
const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m];
return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
}
// Computes rho(r)
inline double RhoOfR(const EAMTableIndex& index, int itype, int jtype) const {
const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m];
return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
}
// Computes the derivative of F(rho)
inline double FPrimeOfRho(const EAMTableIndex& index, int itype) const {
const double* coeff = frho_spline[type2frho[itype]][index.m];
return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
}
// Computes F(rho)
inline double FofRho(const EAMTableIndex& index, int itype) const {
const double* coeff = frho_spline[type2frho[itype]][index.m];
return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
}
// Computes the derivative of z2(r)
inline double Z2PrimeOfR(const EAMTableIndex& index, int itype, int jtype) const {
const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
}
// Computes z2(r)
inline double Z2OfR(const EAMTableIndex& index, int itype, int jtype) const {
const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
}
// Computes pair potential V_ij(r).
inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR) const {
// phi = pair potential energy
// z2 = phi * r
const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
return z2 * oneOverR;
}
// Computes pair potential V_ij(r) and its derivative.
inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR, double& phid) const {
// phi = pair potential energy
// phip = phi'
// z2 = phi * r
// z2p = (phi * r)' = (phi' r) + phi
const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
const double z2p = (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
const double phi = z2 * oneOverR;
phid = z2p * oneOverR - phi * oneOverR;
return phi;
}
// Parameters
// h() polynomial function coefficients
double* hcoeff;
// The number of coefficients in the polynomial.
int nhcoeff;
// This specifies the calculation stage to let the pack/unpack communication routines know
// which data to exchange.
int communicationStage;
};
/// The one-site concentration formulation of CD-EAM.
class PairCDEAM_OneSite : public PairCDEAM
{
public:
/// Constructor.
PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {}
};
/// The two-site concentration formulation of CD-EAM.
class PairCDEAM_TwoSite : public PairCDEAM
{
public:
/// Constructor.
PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {}
};
};
#endif
#endif