lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2011-08-25 15:30:01 +00:00
parent b19a38a92d
commit 9eb5aebed9
7 changed files with 10 additions and 652 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/DSMC/Install.sh
View File

@ -1,15 +0,0 @@
# Install/unInstall package files in LAMMPS
if (test $1 = 1) then
cp pair_dsmc.cpp ..
cp pair_dsmc.h ..
elif (test $1 = 0) then
rm -f ../pair_dsmc.cpp
rm -f ../pair_dsmc.h
fi

525
src/DSMC/pair_dsmc.cpp
View File

@ -1,525 +0,0 @@
/* ----------------------------------------------------------------------
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: Paul Crozier (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_dsmc.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include "update.h"
#include "random_mars.h"
#include "limits.h"
using namespace LAMMPS_NS;
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
/* ---------------------------------------------------------------------- */
PairDSMC::PairDSMC(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
total_number_of_collisions = 0;
max_particles = max_particle_list = 0;
next_particle = NULL;
random = NULL;
}
/* ---------------------------------------------------------------------- */
PairDSMC::~PairDSMC()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(sigma);
memory->destroy(cut);
memory->destroy(V_sigma_max);
memory->destroy(particle_list);
memory->destroy(first);
memory->destroy(number);
}
delete [] next_particle;
delete random;
}
/* ---------------------------------------------------------------------- */
void PairDSMC::compute(int eflag, int vflag)
{
double **x = atom->x;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
for (int i = 1; i <= atom->ntypes; ++i)
for (int j = 0; j < total_ncells; ++j) {
first[i][j] = -1;
number[i][j] = 0;
}
if (atom->nmax > max_particles) {
delete [] next_particle;
max_particles = atom->nmax;
next_particle = new int[max_particles];
}
// find each particle's cell and sort by type
// assume a constant volume and shape simulation domain
// skip particle if outside processor domain
for (int i = 0; i < nlocal; ++i) {
int xcell = static_cast<int>((x[i][0] - domain->boxlo[0])/cellx);
int ycell = static_cast<int>((x[i][1] - domain->boxlo[1])/celly);
int zcell = static_cast<int>((x[i][2] - domain->boxlo[2])/cellz);
if ((xcell < 0) or (xcell > ncellsx-1) or
(ycell < 0) or (ycell > ncellsy-1) or
(zcell < 0) or (zcell > ncellsz-1)) continue;
int icell = xcell + ycell*ncellsx + zcell*ncellsx*ncellsy;
itype = type[i];
next_particle[i] = first[itype][icell];
first[itype][icell] = i;
number[itype][icell]++;
}
for (int icell = 0; icell < total_ncells; ++icell) {
for (itype = 1; itype <= atom->ntypes; ++itype) {
number_of_A = number[itype][icell];
if (number_of_A > max_particle_list) {
max_particle_list = number_of_A;
memory->grow(particle_list,atom->ntypes+1,max_particle_list,
"pair:particle_list");
}
int m = first[itype][icell];
for (int k = 0; k < number_of_A; k++) {
particle_list[itype][k] = m;
m = next_particle[m];
}
}
for (itype = 1; itype <= atom->ntypes; ++itype) {
imass = mass[itype];
number_of_A = number[itype][icell];
for (jtype = itype; jtype <= atom->ntypes; ++jtype) {
jmass = mass[jtype];
number_of_B = number[jtype][icell];
reduced_mass = imass*jmass/(imass + jmass);
total_mass = imass + jmass;
jmass_tmass = jmass/total_mass;
imass_tmass = imass/total_mass;
// if necessary, recompute V_sigma_max values
if (recompute_vsigmamax_stride &&
(update->ntimestep % recompute_vsigmamax_stride == 0))
recompute_V_sigma_max(icell);
// # of collisions to perform for itype-jtype pairs
double &Vs_max = V_sigma_max[itype][jtype];
double num_of_collisions_double = number_of_A * number_of_B *
weighting * Vs_max * update->dt / vol;
if ((itype == jtype) and number_of_B)
num_of_collisions_double *=
0.5 * double(number_of_B - 1) / double(number_of_B);
int num_of_collisions =
convert_double_to_equivalent_int(num_of_collisions_double);
if (num_of_collisions > number_of_A)
error->warning("Pair dsmc: num_of_collisions > number_of_A",0);
if (num_of_collisions > number_of_B)
error->warning("Pair dsmc: num_of_collisions > number_of_B",0);
// perform collisions on pairs of particles in icell
for (int k = 0; k < num_of_collisions; k++) {
if ((number_of_A < 1) or (number_of_B < 1)) break;
int ith_A = static_cast<int>(random->uniform()*number_of_A);
int jth_B = static_cast<int>(random->uniform()*number_of_B);
int i = particle_list[itype][ith_A];
int j = particle_list[jtype][jth_B];
if (i == j) {
k--;
continue;
}
double probability = V_sigma(i,j)/Vs_max;
if (probability > random->uniform()) scatter_random(i,j,icell);
}
}
}
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairDSMC::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(V_sigma_max,n+1,n+1,"pair:V_sigma_max");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairDSMC::settings(int narg, char **arg)
{
if (narg != 6) error->all("Illegal pair_style command");
cut_global = 0.0;
max_cell_size = force->numeric(arg[0]);
seed = force->inumeric(arg[1]);
weighting = force->numeric(arg[2]);
T_ref = force->numeric(arg[3]);
recompute_vsigmamax_stride = force->inumeric(arg[4]);
vsigmamax_samples = force->inumeric(arg[5]);
// initialize Marsaglia RNG with processor-unique seed
if (max_cell_size <= 0.0) error->all("Illegal pair_style command");
if (seed <= 0) error->all("Illegal pair_style command");
if (random) delete random;
random = new RanMars(lmp,seed + comm->me);
kT_ref = force->boltz*T_ref;
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairDSMC::coeff(int narg, char **arg)
{
if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double sigma_one = force->numeric(arg[2]);
double cut_one = cut_global;
if (narg == 4) cut_one = force->numeric(arg[3]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all("Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairDSMC::init_style()
{
ncellsx = ncellsy = ncellsz = 1;
while (((domain->boxhi[0] - domain->boxlo[0])/ncellsx) > max_cell_size)
ncellsx++;
while (((domain->boxhi[1] - domain->boxlo[1])/ncellsy) > max_cell_size)
ncellsy++;
while (((domain->boxhi[2] - domain->boxlo[2])/ncellsz) > max_cell_size)
ncellsz++;
cellx = (domain->boxhi[0] - domain->boxlo[0])/ncellsx;
celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy;
cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz;
if (comm->me == 0) {
if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n",
cellx,celly,cellz);
if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n",
cellx,celly,cellz);
}
total_ncells = ncellsx*ncellsy*ncellsz;
vol = cellx*celly*cellz;
memory->create(particle_list,atom->ntypes+1,0,"pair:particle_list");
memory->create(first,atom->ntypes+1,total_ncells,"pair:first");
memory->create(number,atom->ntypes+1,total_ncells,"pair:number");
for (int i = 1; i <= atom->ntypes; i++)
for (int j = 1; j <= atom->ntypes; j++)
V_sigma_max[i][j] = 0.0;
two_pi = 8.0*atan(1.0);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairDSMC::init_one(int i, int j)
{
if (setflag[i][j] == 0) cut[i][j] = 0.0;
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairDSMC::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairDSMC::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairDSMC::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&max_cell_size,sizeof(double),1,fp);
fwrite(&seed,sizeof(int),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairDSMC::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&max_cell_size,sizeof(double),1,fp);
fread(&seed,sizeof(int),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&max_cell_size,1,MPI_DOUBLE,0,world);
MPI_Bcast(&seed,1,MPI_INT,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
// initialize Marsaglia RNG with processor-unique seed
// same seed that pair_style command initially specified
if (random) delete random;
random = new RanMars(lmp,seed + comm->me);
}
/*-------------------------------------------------------------------------
rezero and recompute the V_sigma_max values this timestep for use during
the next nrezero timesteps
-------------------------------------------------------------------------*/
void PairDSMC::recompute_V_sigma_max(int icell)
{
int i,j,k;
double Vsigma_max = 0;
if (number_of_A && number_of_B) {
for (k = 0; k < vsigmamax_samples; k++) {
i = particle_list[itype]
[static_cast<int>(random->uniform()*number_of_A)];
j = particle_list[jtype]
[static_cast<int>(random->uniform()*number_of_B)];
if (i == j) continue;
Vsigma_max = MAX(Vsigma_max,V_sigma(i,j));
}
}
V_sigma_max[itype][jtype] = Vsigma_max;
}
/*-------------------------------------------------------------------------
VHS model
compute the velocity vector difference between i and j and multiply by
their combined collision cross section, sigma, for neutral-neutral
collisions using the Variable Hard Sphere model
-------------------------------------------------------------------------*/
double PairDSMC::V_sigma(int i, int j)
{
double relative_velocity_sq,relative_velocity,pair_sigma;
double delv[3];
double *vi = atom->v[i];
double *vj = atom->v[j];
subtract3d(vi,vj,delv);
relative_velocity_sq = dot3d(delv,delv);
relative_velocity = sqrt(relative_velocity_sq);
// from Bird eq 4.63, and omega=0.67
// (omega - 0.5) = 0.17
// 1/GAMMA(2.5 - omega) = 1.06418029298371
if (relative_velocity_sq != 0.0)
pair_sigma = sigma[itype][jtype]*
pow(kT_ref/(0.5*reduced_mass*relative_velocity_sq),0.17) *
1.06418029298371;
else
pair_sigma = 0.0;
return relative_velocity*pair_sigma;
}
/*-------------------------------------------------------------------------
generate new velocities for collided particles
-------------------------------------------------------------------------*/
void PairDSMC::scatter_random(int i, int j, int icell)
{
double mag_delv,cos_phi,cos_squared,r,theta;
double delv[3],vcm[3];
double *vi = atom->v[i];
double *vj = atom->v[j];
subtract3d(vi,vj,delv);
if (itype == jtype) mag_delv = sqrt(dot3d(delv,delv))*0.5;
else mag_delv = sqrt(dot3d(delv,delv));
cos_phi = 1.0 - (2.0*random->uniform());
cos_squared = MIN(1.0,cos_phi*cos_phi);
r = sqrt(1.0 - cos_squared);
delv[0] = cos_phi*mag_delv;
theta = two_pi*random->uniform();
delv[1] = r*mag_delv*cos(theta);
delv[2] = r*mag_delv*sin(theta);
if (itype == jtype) {
vcm[0] = (vi[0]+vj[0])*0.5;
vcm[1] = (vi[1]+vj[1])*0.5;
vcm[2] = (vi[2]+vj[2])*0.5;
vi[0] = vcm[0] + delv[0];
vi[1] = vcm[1] + delv[1];
vi[2] = vcm[2] + delv[2];
vj[0] = vcm[0] - delv[0];
vj[1] = vcm[1] - delv[1];
vj[2] = vcm[2] - delv[2];
} else {
vcm[0] = vi[0]*imass_tmass + vj[0]*jmass_tmass;
vcm[1] = vi[1]*imass_tmass + vj[1]*jmass_tmass;
vcm[2] = vi[2]*imass_tmass + vj[2]*jmass_tmass;
vi[0] = vcm[0] + delv[0]*jmass_tmass;
vi[1] = vcm[1] + delv[1]*jmass_tmass;
vi[2] = vcm[2] + delv[2]*jmass_tmass;
vj[0] = vcm[0] - delv[0]*imass_tmass;
vj[1] = vcm[1] - delv[1]*imass_tmass;
vj[2] = vcm[2] - delv[2]*imass_tmass;
}
total_number_of_collisions++;
}
/* ----------------------------------------------------------------------
This method converts the double supplied by the calling function into
an int, which is returned. By adding a random number between 0 and 1
to the double before converting it to an int, we ensure that,
statistically, we round down with probability identical to the
remainder and up the rest of the time. So even though we're using an
integer, we're statistically matching the exact expression represented
by the double.
------------------------------------------------------------------------- */
int PairDSMC::convert_double_to_equivalent_int(double input_double)
{
if (input_double > INT_MAX)
error->all("Tried to convert a double to int, but input_double > INT_MAX");
int output_int = static_cast<int>(input_double + random->uniform());
return output_int;
}

109
src/DSMC/pair_dsmc.h
View File

@ -1,109 +0,0 @@
/* ----------------------------------------------------------------------
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(dsmc,PairDSMC)
#else
#ifndef LMP_PAIR_DSMC_H
#define LMP_PAIR_DSMC_H
#include "pair.h"
namespace LAMMPS_NS {
class PairDSMC : public Pair {
public:
PairDSMC(class LAMMPS *);
virtual ~PairDSMC();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
private:
double cut_global;
double **cut;
double **sigma;
double cellx;
double celly;
double cellz;
int ncellsx;
int ncellsy;
int ncellsz;
int total_ncells;
int total_number_of_collisions;
int recompute_vsigmamax_stride;
int vsigmamax_samples;
double T_ref;
double kT_ref;
double two_pi;
double max_cell_size;
int seed;
int number_of_A;
int number_of_B;
int max_particle_list;
class RanMars *random;
int **particle_list;
int **first;
int **number;
double **V_sigma_max;
int max_particles;
int *next_particle;
int itype;
int jtype;
double imass;
double jmass;
double total_mass;
double reduced_mass;
double imass_tmass;
double jmass_tmass;
double vol;
double weighting;
void allocate();
void recompute_V_sigma_max(int);
double V_sigma(int, int);
void scatter_random(int, int, int);
int convert_double_to_equivalent_int(double);
inline void subtract3d(const double *v1, const double *v2, double *v3) {
v3[0] = v2[0] - v1[0];
v3[1] = v2[1] - v1[1];
v3[2] = v2[2] - v1[2];
}
inline double dot3d(const double *v1, const double *v2) {
return( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] );
}
};
}
#endif
#endif

4
src/Makefile
View File

@ -13,8 +13,8 @@ OBJ = $(SRC:.cpp=.o)
# Package variables
PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
kspace manybody meam molecule opt peri poems reax replica \
PACKAGE = asphere class2 colloid dipole gpu granular \
kspace manybody mc meam molecule opt peri poems reax replica \
shock srd xtc
PACKUSER = user-misc user-atc user-awpmd user-cg-cmm \

2
src/fix_restrain.cpp
View File

@ -288,7 +288,7 @@ void FixRestrain::restrain_dihedral()
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Restrain problem: %d %d %d %d %d %d",
sprintf(str,"Restrain problem: %d " BIGINT_FORMAT " %d %d %d %d",
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(str);

1
src/force.h
View File

@ -21,6 +21,7 @@ namespace LAMMPS_NS {
class Force : protected Pointers {
public:
double boltz; // Boltzmann constant (eng/degree-K)
double hplanck; // Planck's constant (energy-time)
double mvv2e; // conversion of mv^2 to energy
double ftm2v; // conversion of ft/m to velocity
double mv2d; // conversion of mass/volume to density

6
src/update.cpp
View File

@ -118,6 +118,7 @@ void Update::set_units(const char *style)
if (strcmp(style,"lj") == 0) {
force->boltz = 1.0;
force->hplanck = 0.18292026; // using LJ parameters for argon
force->mvv2e = 1.0;
force->ftm2v = 1.0;
force->mv2d = 1.0;
@ -135,6 +136,7 @@ void Update::set_units(const char *style)
} else if (strcmp(style,"real") == 0) {
force->boltz = 0.0019872067;
force->hplanck = 95.306976368;
force->mvv2e = 48.88821291 * 48.88821291;
force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
force->mv2d = 1.0 / 0.602214179;
@ -152,6 +154,7 @@ void Update::set_units(const char *style)
} else if (strcmp(style,"metal") == 0) {
force->boltz = 8.617343e-5;
force->hplanck = 4.135667403e-3;
force->mvv2e = 1.0364269e-4;
force->ftm2v = 1.0 / 1.0364269e-4;
force->mv2d = 1.0 / 0.602214179;
@ -169,6 +172,7 @@ void Update::set_units(const char *style)
} else if (strcmp(style,"si") == 0) {
force->boltz = 1.3806504e-23;
force->hplanck = 6.62606896e-34;
force->mvv2e = 1.0;
force->ftm2v = 1.0;
force->mv2d = 1.0;
@ -186,6 +190,7 @@ void Update::set_units(const char *style)
} else if (strcmp(style,"cgs") == 0) {
force->boltz = 1.3806504e-16;
force->hplanck = 6.62606896e-27;
force->mvv2e = 1.0;
force->ftm2v = 1.0;
force->mv2d = 1.0;
@ -203,6 +208,7 @@ void Update::set_units(const char *style)
} else if (strcmp(style,"electron") == 0) {
force->boltz = 3.16681534e-6;
force->hplanck = 0.1519829846;
force->mvv2e = 1.06657236;
force->ftm2v = 0.937582899;
force->mv2d = 1.0;