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@10200 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-06-29 21:28:36 +00:00
parent b0007d45b7
commit d2c3a96a9a
8 changed files with 1250 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

323
src/USER-OMP/pair_eam_alloy_omp.cpp.orig Normal file
View File

@ -0,0 +1,323 @@
/* ----------------------------------------------------------------------
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 authors: Stephen Foiles (SNL), Murray Daw (SNL)
------------------------------------------------------------------------- */
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_eam_alloy_omp.h"
#include "atom.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
PairEAMAlloyOMP::PairEAMAlloyOMP(LAMMPS *lmp) : PairEAMOMP(lmp)
{
one_coeff = 1;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
read DYNAMO setfl file
------------------------------------------------------------------------- */
void PairEAMAlloyOMP::coeff(int narg, char **arg)
{
int i,j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read EAM setfl file
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
delete [] setfl->mass;
memory->destroy(setfl->frho);
memory->destroy(setfl->rhor);
memory->destroy(setfl->z2r);
delete setfl;
}
setfl = new Setfl();
read_file(arg[2]);
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < setfl->nelements; j++)
if (strcmp(arg[i],setfl->elements[j]) == 0) break;
if (j < setfl->nelements) map[i-2] = j;
else error->all(FLERR,"No matching element in EAM potential file");
}
// clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes;
for (i = 1; i <= n; i++)
for (j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass of atom type if i = j
int count = 0;
for (i = 1; i <= n; i++) {
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(i,setfl->mass[map[i]]);
count++;
}
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
read a multi-element DYNAMO setfl file
------------------------------------------------------------------------- */
void PairEAMAlloyOMP::read_file(char *filename)
{
Setfl *file = setfl;
// open potential file
int me = comm->me;
FILE *fptr;
char line[MAXLINE];
if (me == 0) {
fptr = fopen(filename,"r");
if (fptr == NULL) {
char str[128];
sprintf(str,"Cannot open EAM potential file %s",filename);
error->one(FLERR,str);
}
}
// read and broadcast header
// extract element names from nelements line
int n;
if (me == 0) {
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
sscanf(line,"%d",&file->nelements);
int nwords = atom->count_words(line);
if (nwords != file->nelements + 1)
error->all(FLERR,"Incorrect element names in EAM potential file");
char **words = new char*[file->nelements+1];
nwords = 0;
strtok(line," \t\n\r\f");
while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
file->elements = new char*[file->nelements];
for (int i = 0; i < file->nelements; i++) {
n = strlen(words[i]) + 1;
file->elements[i] = new char[n];
strcpy(file->elements[i],words[i]);
}
delete [] words;
if (me == 0) {
fgets(line,MAXLINE,fptr);
sscanf(line,"%d %lg %d %lg %lg",
&file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
}
MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
MPI_Bcast(&file->nr,1,MPI_INT,0,world);
MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
file->mass = new double[file->nelements];
memory->create(file->frho,file->nelements,file->nrho+1,"pair:frho");
memory->create(file->rhor,file->nelements,file->nr+1,"pair:rhor");
memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,
"pair:z2r");
int i,j,tmp;
for (i = 0; i < file->nelements; i++) {
if (me == 0) {
fgets(line,MAXLINE,fptr);
sscanf(line,"%d %lg",&tmp,&file->mass[i]);
}
MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
if (me == 0) grab(fptr,file->nr,&file->rhor[i][1]);
MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world);
}
for (i = 0; i < file->nelements; i++)
for (j = 0; j <= i; j++) {
if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
}
// close the potential file
if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
copy read-in setfl potential to standard array format
------------------------------------------------------------------------- */
void PairEAMAlloyOMP::file2array()
{
int i,j,m,n;
int ntypes = atom->ntypes;
// set function params directly from setfl file
nrho = setfl->nrho;
nr = setfl->nr;
drho = setfl->drho;
dr = setfl->dr;
// ------------------------------------------------------------------
// setup frho arrays
// ------------------------------------------------------------------
// allocate frho arrays
// nfrho = # of setfl elements + 1 for zero array
nfrho = setfl->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
// copy each element's frho to global frho
for (i = 0; i < setfl->nelements; i++)
for (m = 1; m <= nrho; m++) frho[i][m] = setfl->frho[i][m];
// add extra frho of zeroes for non-EAM types to point to (pair hybrid)
// this is necessary b/c fp is still computed for non-EAM atoms
for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0;
// type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to
// if atom type doesn't point to element (non-EAM atom in pair hybrid)
// then map it to last frho array of zeroes
for (i = 1; i <= ntypes; i++)
if (map[i] >= 0) type2frho[i] = map[i];
else type2frho[i] = nfrho-1;
// ------------------------------------------------------------------
// setup rhor arrays
// ------------------------------------------------------------------
// allocate rhor arrays
// nrhor = # of setfl elements
nrhor = setfl->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
// copy each element's rhor to global rhor
for (i = 0; i < setfl->nelements; i++)
for (m = 1; m <= nr; m++) rhor[i][m] = setfl->rhor[i][m];
// type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to
// for setfl files, I,J mapping only depends on I
// OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++)
type2rhor[i][j] = map[i];
// ------------------------------------------------------------------
// setup z2r arrays
// ------------------------------------------------------------------
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of setfl elements
nz2r = setfl->nelements * (setfl->nelements+1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
n = 0;
for (i = 0; i < setfl->nelements; i++)
for (j = 0; j <= i; j++) {
for (m = 1; m <= nr; m++) z2r[n][m] = setfl->z2r[i][j][m];
n++;
}
// type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to
// set of z2r arrays only fill lower triangular Nelement matrix
// value = n = sum over rows of lower-triangular matrix until reach irow,icol
// swap indices when irow < icol to stay lower triangular
// if map = -1 (non-EAM atom in pair hybrid):
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];
icol = map[j];
if (irow == -1 || icol == -1) {
type2z2r[i][j] = 0;
continue;
}
if (irow < icol) {
irow = map[j];
icol = map[i];
}
n = 0;
for (m = 0; m < irow; m++) n += m + 1;
n += icol;
type2z2r[i][j] = n;
}
}
}

332
src/USER-OMP/pair_eam_fs_omp.cpp.orig Normal file
View File

@ -0,0 +1,332 @@
/* ----------------------------------------------------------------------
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 authors: Tim Lau (MIT)
------------------------------------------------------------------------- */
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_eam_fs_omp.h"
#include "atom.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
PairEAMFSOMP::PairEAMFSOMP(LAMMPS *lmp) : PairEAMOMP(lmp)
{
one_coeff = 1;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
read EAM Finnis-Sinclair file
------------------------------------------------------------------------- */
void PairEAMFSOMP::coeff(int narg, char **arg)
{
int i,j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read EAM Finnis-Sinclair file
if (fs) {
for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i];
delete [] fs->elements;
delete [] fs->mass;
memory->destroy(fs->frho);
memory->destroy(fs->rhor);
memory->destroy(fs->z2r);
delete fs;
}
fs = new Fs();
read_file(arg[2]);
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < fs->nelements; j++)
if (strcmp(arg[i],fs->elements[j]) == 0) break;
if (j < fs->nelements) map[i-2] = j;
else error->all(FLERR,"No matching element in EAM potential file");
}
// clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes;
for (i = 1; i <= n; i++)
for (j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass of atom type if i = j
int count = 0;
for (i = 1; i <= n; i++) {
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(i,fs->mass[map[i]]);
count++;
}
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
read a multi-element DYNAMO setfl file
------------------------------------------------------------------------- */
void PairEAMFSOMP::read_file(char *filename)
{
Fs *file = fs;
// open potential file
int me = comm->me;
FILE *fptr;
char line[MAXLINE];
if (me == 0) {
fptr = fopen(filename,"r");
if (fptr == NULL) {
char str[128];
sprintf(str,"Cannot open EAM potential file %s",filename);
error->one(FLERR,str);
}
}
// read and broadcast header
// extract element names from nelements line
int n;
if (me == 0) {
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
fgets(line,MAXLINE,fptr);
n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
sscanf(line,"%d",&file->nelements);
int nwords = atom->count_words(line);
if (nwords != file->nelements + 1)
error->all(FLERR,"Incorrect element names in EAM potential file");
char **words = new char*[file->nelements+1];
nwords = 0;
strtok(line," \t\n\r\f");
while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
file->elements = new char*[file->nelements];
for (int i = 0; i < file->nelements; i++) {
n = strlen(words[i]) + 1;
file->elements[i] = new char[n];
strcpy(file->elements[i],words[i]);
}
delete [] words;
if (me == 0) {
fgets(line,MAXLINE,fptr);
sscanf(line,"%d %lg %d %lg %lg",
&file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
}
MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
MPI_Bcast(&file->nr,1,MPI_INT,0,world);
MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
file->mass = new double[file->nelements];
memory->create(file->frho,file->nelements,file->nrho+1,
"pair:frho");
memory->create(file->rhor,file->nelements,file->nelements,
file->nr+1,"pair:rhor");
memory->create(file->z2r,file->nelements,file->nelements,
file->nr+1,"pair:z2r");
int i,j,tmp;
for (i = 0; i < file->nelements; i++) {
if (me == 0) {
fgets(line,MAXLINE,fptr);
sscanf(line,"%d %lg",&tmp,&file->mass[i]);
}
MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
for (j = 0; j < file->nelements; j++) {
if (me == 0) grab(fptr,file->nr,&file->rhor[i][j][1]);
MPI_Bcast(&file->rhor[i][j][1],file->nr,MPI_DOUBLE,0,world);
}
}
for (i = 0; i < file->nelements; i++)
for (j = 0; j <= i; j++) {
if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
}
// close the potential file
if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
copy read-in setfl potential to standard array format
------------------------------------------------------------------------- */
void PairEAMFSOMP::file2array()
{
int i,j,m,n;
int ntypes = atom->ntypes;
// set function params directly from fs file
nrho = fs->nrho;
nr = fs->nr;
drho = fs->drho;
dr = fs->dr;
// ------------------------------------------------------------------
// setup frho arrays
// ------------------------------------------------------------------
// allocate frho arrays
// nfrho = # of fs elements + 1 for zero array
nfrho = fs->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
// copy each element's frho to global frho
for (i = 0; i < fs->nelements; i++)
for (m = 1; m <= nrho; m++) frho[i][m] = fs->frho[i][m];
// add extra frho of zeroes for non-EAM types to point to (pair hybrid)
// this is necessary b/c fp is still computed for non-EAM atoms
for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0;
// type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to
// if atom type doesn't point to element (non-EAM atom in pair hybrid)
// then map it to last frho array of zeroes
for (i = 1; i <= ntypes; i++)
if (map[i] >= 0) type2frho[i] = map[i];
else type2frho[i] = nfrho-1;
// ------------------------------------------------------------------
// setup rhor arrays
// ------------------------------------------------------------------
// allocate rhor arrays
// nrhor = square of # of fs elements
nrhor = fs->nelements * fs->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
// copy each element pair rhor to global rhor
n = 0;
for (i = 0; i < fs->nelements; i++)
for (j = 0; j < fs->nelements; j++) {
for (m = 1; m <= nr; m++) rhor[n][m] = fs->rhor[i][j][m];
n++;
}
// type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to
// for fs files, there is a full NxN set of rhor arrays
// OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++)
type2rhor[i][j] = map[i] * fs->nelements + map[j];
// ------------------------------------------------------------------
// setup z2r arrays
// ------------------------------------------------------------------
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of fs elements
nz2r = fs->nelements * (fs->nelements+1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
n = 0;
for (i = 0; i < fs->nelements; i++)
for (j = 0; j <= i; j++) {
for (m = 1; m <= nr; m++) z2r[n][m] = fs->z2r[i][j][m];
n++;
}
// type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to
// set of z2r arrays only fill lower triangular Nelement matrix
// value = n = sum over rows of lower-triangular matrix until reach irow,icol
// swap indices when irow < icol to stay lower triangular
// if map = -1 (non-EAM atom in pair hybrid):
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];
icol = map[j];
if (irow == -1 || icol == -1) {
type2z2r[i][j] = 0;
continue;
}
if (irow < icol) {
irow = map[j];
icol = map[i];
}
n = 0;
for (m = 0; m < irow; m++) n += m + 1;
n += icol;
type2z2r[i][j] = n;
}
}
}

22
src/USER-OMP/pair_gran_hertz_history_omp.cpp
View File

@ -18,6 +18,7 @@
#include "comm.h" #include "comm.h"
#include "fix.h" #include "fix.h"
#include "force.h" #include "force.h"
#include "memory.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "update.h" #include "update.h"
@ -49,12 +50,23 @@ void PairGranHertzHistoryOMP::compute(int eflag, int vflag)
computeflag = 1; computeflag = 1;
const int shearupdate = (update->setupflag) ? 0 : 1; const int shearupdate = (update->setupflag) ? 0 : 1;
// update body ptr and values for ghost atoms if using FixRigid masses // update rigid body info for owned & ghost atoms if using FixRigid masses
// body[i] = which body atom I is in, -1 if none
// mass_body = mass of each rigid body
if (fix_rigid && neighbor->ago == 0) { if (fix_rigid && neighbor->ago == 0) {
int tmp; int tmp;
body = (int *) fix_rigid->extract("body",tmp); int *body = (int *) fix_rigid->extract("body",tmp);
mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
if (atom->nmax > nmax) {
memory->destroy(mass_rigid);
nmax = atom->nmax;
memory->create(mass_rigid,nmax,"pair:mass_rigid");
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
else mass_rigid[i] = 0.0;
comm->forward_comm_pair(this); comm->forward_comm_pair(this);
} }
@ -193,8 +205,8 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
mj = mass[type[j]]; mj = mass[type[j]];
} }
if (fix_rigid) { if (fix_rigid) {
if (body[i] >= 0) mi = mass_rigid[body[i]]; if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
if (body[j] >= 0) mj = mass_rigid[body[j]]; if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
} }
meff = mi*mj / (mi+mj); meff = mi*mj / (mi+mj);

22
src/USER-OMP/pair_gran_hooke_history_omp.cpp
View File

@ -18,6 +18,7 @@
#include "comm.h" #include "comm.h"
#include "fix.h" #include "fix.h"
#include "force.h" #include "force.h"
#include "memory.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "update.h" #include "update.h"
@ -50,12 +51,23 @@ void PairGranHookeHistoryOMP::compute(int eflag, int vflag)
computeflag = 1; computeflag = 1;
const int shearupdate = (update->setupflag) ? 0 : 1; const int shearupdate = (update->setupflag) ? 0 : 1;
// update body ptr and values for ghost atoms if using FixRigid masses // update rigid body info for owned & ghost atoms if using FixRigid masses
// body[i] = which body atom I is in, -1 if none
// mass_body = mass of each rigid body
if (fix_rigid && neighbor->ago == 0) { if (fix_rigid && neighbor->ago == 0) {
int tmp; int tmp;
body = (int *) fix_rigid->extract("body",tmp); int *body = (int *) fix_rigid->extract("body",tmp);
mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
if (atom->nmax > nmax) {
memory->destroy(mass_rigid);
nmax = atom->nmax;
memory->create(mass_rigid,nmax,"pair:mass_rigid");
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
else mass_rigid[i] = 0.0;
comm->forward_comm_pair(this); comm->forward_comm_pair(this);
} }
@ -198,8 +210,8 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
mj = mass[type[j]]; mj = mass[type[j]];
} }
if (fix_rigid) { if (fix_rigid) {
if (body[i] >= 0) mi = mass_rigid[body[i]]; if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
if (body[j] >= 0) mj = mass_rigid[body[j]]; if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
} }
meff = mi*mj / (mi+mj); meff = mi*mj / (mi+mj);

22
src/USER-OMP/pair_gran_hooke_omp.cpp
View File

@ -18,6 +18,7 @@
#include "comm.h" #include "comm.h"
#include "fix.h" #include "fix.h"
#include "force.h" #include "force.h"
#include "memory.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
@ -45,12 +46,23 @@ void PairGranHookeOMP::compute(int eflag, int vflag)
const int nthreads = comm->nthreads; const int nthreads = comm->nthreads;
const int inum = list->inum; const int inum = list->inum;
// update body ptr and values for ghost atoms if using FixRigid masses // update rigid body info for owned & ghost atoms if using FixRigid masses
// body[i] = which body atom I is in, -1 if none
// mass_body = mass of each rigid body
if (fix_rigid && neighbor->ago == 0) { if (fix_rigid && neighbor->ago == 0) {
int tmp; int tmp;
body = (int *) fix_rigid->extract("body",tmp); int *body = (int *) fix_rigid->extract("body",tmp);
mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
if (atom->nmax > nmax) {
memory->destroy(mass_rigid);
nmax = atom->nmax;
memory->create(mass_rigid,nmax,"pair:mass_rigid");
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
else mass_rigid[i] = 0.0;
comm->forward_comm_pair(this); comm->forward_comm_pair(this);
} }
@ -172,8 +184,8 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr)
mj = mass[type[j]]; mj = mass[type[j]];
} }
if (fix_rigid) { if (fix_rigid) {
if (body[i] >= 0) mi = mass_rigid[body[i]]; if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
if (body[j] >= 0) mj = mass_rigid[body[j]]; if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
} }
meff = mi*mj / (mi+mj); meff = mi*mj / (mi+mj);

199
src/USER-OMP/pair_nb3b_harmonic_omp.cpp Normal file
View File

@ -0,0 +1,199 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "math.h"
#include "pair_nb3b_harmonic_omp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "suffix.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairNb3bHarmonicOMP::PairNb3bHarmonicOMP(LAMMPS *lmp) :
PairNb3bHarmonic(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairNb3bHarmonicOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
if (evflag) {
if (eflag) {
eval<1,1>(ifrom, ito, thr);
} else {
eval<1,0>(ifrom, ito, thr);
}
} else eval<0,0>(ifrom, ito, thr);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG>
void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,kk,jnum,jnumm1,itag,jtag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fj[3],fk[3];
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
const int * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double fxtmp,fytmp,fztmp;
// loop over full neighbor list of my atoms
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i].x;
ytmp = x[i].y;
ztmp = x[i].z;
fxtmp = fytmp = fztmp = 0.0;
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j].z < ztmp) continue;
if (x[j].z == ztmp && x[j].y < ytmp) continue;
if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
}
jtype = map[type[j]];
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx*delx + dely*dely + delz*delz;
ijparam = elem2param[itype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue;
}
jnumm1 = jnum - 1;
for (jj = 0; jj < jnumm1; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
delr1[0] = x[j].x - xtmp;
delr1[1] = x[j].y - ytmp;
delr1[2] = x[j].z - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > params[ijparam].cutsq) continue;
double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0;
for (kk = jj+1; kk < jnum; kk++) {
k = jlist[kk];
k &= NEIGHMASK;
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
delr2[0] = x[k].x - xtmp;
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > params[ikparam].cutsq) continue;
threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
fxtmp -= fj[0] + fk[0];
fytmp -= fj[1] + fk[1];
fztmp -= fj[2] + fk[2];
fjxtmp += fj[0];
fjytmp += fj[1];
fjztmp += fj[2];
f[k].x += fk[0];
f[k].y += fk[1];
f[k].z += fk[2];
if (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;
f[j].y += fjytmp;
f[j].z += fjztmp;
}
f[i].x += fxtmp;
f[i].y += fytmp;
f[i].z += fztmp;
}
}
/* ---------------------------------------------------------------------- */
double PairNb3bHarmonicOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairNb3bHarmonic::memory_usage();
return bytes;
}

48
src/USER-OMP/pair_nb3b_harmonic_omp.h Normal file
View File

@ -0,0 +1,48 @@
/* -*- 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: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(nb3b/harmonic/omp,PairNb3bHarmonicOMP)
#else
#ifndef LMP_PAIR_NB3BHARMONIC_OMP_H
#define LMP_PAIR_NB3BHARMONIC_OMP_H
#include "pair_nb3b_harmonic.h"
#include "thr_omp.h"
namespace LAMMPS_NS {
class PairNb3bHarmonicOMP : public PairNb3bHarmonic, public ThrOMP {
public:
PairNb3bHarmonicOMP(class LAMMPS *);
virtual void compute(int, int);
virtual double memory_usage();
private:
template <int EVFLAG, int EFLAG>
void eval(int ifrom, int ito, ThrData * const thr);
};
}
#endif
#endif

297
src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig Normal file
View File

@ -0,0 +1,297 @@
/* ----------------------------------------------------------------------
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_omp.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"
#include "math_special.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathSpecial;
#define MAXLINE 1024
#define DELTA 4
/* ----------------------------------------------------------------------
Fermi-like smoothing function
------------------------------------------------------------------------- */
static double F_fermi(const double r, const double expsc, const double cut)
{
return 1.0 / (1.0 + exp(-expsc*(r-cut)));
}
/* ----------------------------------------------------------------------
Fermi-like smoothing function derivative with respect to r
------------------------------------------------------------------------- */
static double F_fermi_d(const double r, const double expsc, const double cut)
{
return expsc*exp(-expsc*(r-cut)) / square(1.0 + exp(-expsc*(r-cut)));
}
/* ---------------------------------------------------------------------- */
PairTersoffZBLOMP::PairTersoffZBLOMP(LAMMPS *lmp) : PairTersoffOMP(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");
}
/* ---------------------------------------------------------------------- */
void PairTersoffZBLOMP::read_file(char *file)
{
int params_per_line = 21;
char **words = new char*[params_per_line+1];
memory->sfree(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 PairTersoffZBLOMP::force_zeta(Param *param, double rsq, double zeta_ij,
double &fforce, double &prefactor,
int eflag, double &eng)
{
double r,fa,fa_d,bij;
r = sqrt(rsq);
fa = (r > param->bigr + param->bigd) ? 0.0 :
-param->bigb * exp(-param->lam2 * r) * ters_fc(r,param) *
F_fermi(r,param->ZBLexpscale,param->ZBLcut);
fa_d = (r > param->bigr + param->bigd) ? 0.0 :
param->bigb * exp(-param->lam2 * r) *
(param->lam2 * ters_fc(r,param) *
F_fermi(r,param->ZBLexpscale,param->ZBLcut) -
ters_fc_d(r,param) * F_fermi(r,param->ZBLexpscale,param->ZBLcut)
- ters_fc(r,param) * F_fermi_d(r,param->ZBLexpscale,param->ZBLcut));
bij = ters_bij(zeta_ij,param);
fforce = 0.5*bij*fa_d / r;
prefactor = -0.5*fa * ters_bij_d(zeta_ij,param);
if (eflag) eng = 0.5*bij*fa;
}
/* ---------------------------------------------------------------------- */
void PairTersoffZBLOMP::repulsive(Param *param, double rsq, double &fforce,
int eflag, double &eng)
{
double r,tmp_fc,tmp_fc_d,tmp_exp;
// Tersoff repulsive portion
r = sqrt(rsq);
tmp_fc = ters_fc(r,param);
tmp_fc_d = ters_fc_d(r,param);
tmp_exp = exp(-param->lam1 * r);
double fforce_ters = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1);
double eng_ters = tmp_fc * param->biga * tmp_exp;
// ZBL repulsive portion
double esq = square(global_e);
double a_ij = (0.8854*global_a_0) /
(pow(param->Z_i,0.23) + pow(param->Z_j,0.23));
double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0);
double r_ov_a = r/a_ij;
double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) +
0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);
double dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) -
0.9423*0.5099*exp(-0.9423*r_ov_a) -
0.4029*0.2802*exp(-0.4029*r_ov_a) -
0.2016*0.02817*exp(-0.2016*r_ov_a));
double fforce_ZBL = premult*-rsq* phi + premult/r*dphi;
double eng_ZBL = premult/r*phi;
// combine two parts with smoothing by Fermi-like function
fforce = -(-F_fermi_d(r,param->ZBLexpscale,param->ZBLcut) * eng_ZBL +
(1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*fforce_ZBL +
F_fermi_d(r,param->ZBLexpscale,param->ZBLcut)*eng_ters +
F_fermi(r,param->ZBLexpscale,param->ZBLcut)*fforce_ters) / r;
if (eflag)
eng = (1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*eng_ZBL +
F_fermi(r,param->ZBLexpscale,param->ZBLcut)*eng_ters;
}