diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp
index 4724d516c..a8138067d 100644
--- a/src/GPU/pppm_gpu.cpp
+++ b/src/GPU/pppm_gpu.cpp
@@ -1,984 +1,714 @@
/* ----------------------------------------------------------------------
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: Mike Brown (ORNL), Axel Kohlmeyer (Temple)
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "pppm_gpu.h"
#include "atom.h"
#include "comm.h"
+#include "commgrid.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "domain.h"
#include "fft3d_wrap.h"
#include "remap_wrap.h"
#include "gpu_extra.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include "universe.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXORDER 7
#define OFFSET 16384
#define SMALL 0.00001
#define LARGE 10000.0
#define EPS_HOC 1.0e-7
+enum{REVERSE_RHO};
+enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
+
#ifdef FFT_SINGLE
#define ZEROF 0.0f
#define ONEF 1.0f
#else
#define ZEROF 0.0
#define ONEF 1.0
#endif
-// External functions from cuda library for atom decomposition
+// external functions from cuda library for atom decomposition
#ifdef FFT_SINGLE
#define PPPM_GPU_API(api) pppm_gpu_ ## api ## _f
#else
#define PPPM_GPU_API(api) pppm_gpu_ ## api ## _d
#endif
FFT_SCALAR* PPPM_GPU_API(init)(const int nlocal, const int nall, FILE *screen,
const int order, const int nxlo_out,
const int nylo_out, const int nzlo_out,
const int nxhi_out, const int nyhi_out,
const int nzhi_out, FFT_SCALAR **rho_coeff,
FFT_SCALAR **_vd_brick,
const double slab_volfactor,
const int nx_pppm, const int ny_pppm,
const int nz_pppm, const bool split,
int &success);
void PPPM_GPU_API(clear)(const double poisson_time);
int PPPM_GPU_API(spread)(const int ago, const int nlocal, const int nall,
double **host_x, int *host_type, bool &success,
double *host_q, double *boxlo, const double delxinv,
const double delyinv, const double delzinv);
void PPPM_GPU_API(interp)(const FFT_SCALAR qqrd2e_scale);
double PPPM_GPU_API(bytes)();
void PPPM_GPU_API(forces)(double **f);
/* ---------------------------------------------------------------------- */
PPPMGPU::PPPMGPU(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg)
{
if (narg != 1) error->all(FLERR,"Illegal kspace_style pppm/gpu command");
density_brick_gpu = vd_brick = NULL;
kspace_split = false;
im_real_space = false;
GPU_EXTRA::gpu_ready(lmp->modify, lmp->error);
}
/* ----------------------------------------------------------------------
free all memory
------------------------------------------------------------------------- */
PPPMGPU::~PPPMGPU()
{
PPPM_GPU_API(clear)(poisson_time);
}
/* ----------------------------------------------------------------------
called once before run
------------------------------------------------------------------------- */
void PPPMGPU::init()
{
+ // PPPM init manages all arrays except density_brick_gpu and vd_brick
+ // thru its deallocate(), allocate()
+ // NOTE: could free density_brick and vdxyz_brick after PPPM allocates them,
+ // before allocating db_gpu and vd_brick down below, if don't need,
+ // if do this, make sure to set them to NULL
+ // NOTE: delete/realloc of cg necessary b/c packing 4 values per grid pt,
+ // not 3 as PPPM does - probably a better way to account for this
+ // in PPPM::init()
+
+ destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
+ destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
+ density_brick_gpu = vd_brick = NULL;
+
PPPM::init();
+ if (differentiation_flag == 0) {
+ delete cg;
+ int (*procneigh)[2] = comm->procneigh;
+ cg = new CommGrid(lmp,world,4,1,
+ nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
+ nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out,
+ procneigh[0][0],procneigh[0][1],procneigh[1][0],
+ procneigh[1][1],procneigh[2][0],procneigh[2][1]);
+ }
+
+ // unsupported option
+
if (differentiation_flag == 1)
- error->all(FLERR,"Cannot (yet) do analytic differentiation with pppm/gpu.");
+ error->all(FLERR,"Cannot (yet) do analytic differentiation with pppm/gpu");
if (strcmp(update->integrate_style,"verlet/split") == 0) {
kspace_split=true;
old_nlocal = 0;
}
if (kspace_split && universe->iworld == 0) {
im_real_space = true;
return;
}
// GPU precision specific init
if (order>8)
error->all(FLERR,"Cannot use order greater than 8 with pppm/gpu.");
PPPM_GPU_API(clear)(poisson_time);
int success;
FFT_SCALAR *data, *h_brick;
h_brick = PPPM_GPU_API(init)(atom->nlocal, atom->nlocal+atom->nghost, screen,
order, nxlo_out, nylo_out, nzlo_out, nxhi_out,
nyhi_out, nzhi_out, rho_coeff, &data,
slab_volfactor,nx_pppm,ny_pppm,nz_pppm,
kspace_split,success);
GPU_EXTRA::check_flag(success,error,world);
+ // allocate density_brick_gpu and vd_brick
+
density_brick_gpu =
create_3d_offset(nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"pppm:density_brick_gpu",h_brick,1);
vd_brick =
create_3d_offset(nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"pppm:vd_brick",data,4);
- poisson_time=0;
+ poisson_time = 0.0;
}
/* ----------------------------------------------------------------------
compute the PPPMGPU long-range force, energy, virial
------------------------------------------------------------------------- */
void PPPMGPU::compute(int eflag, int vflag)
{
int nago;
if (kspace_split) {
if (im_real_space) return;
if (atom->nlocal > old_nlocal) {
nago=0;
old_nlocal = atom->nlocal;
- } else
- nago=1;
- } else
- nago=neighbor->ago;
+ } else nago = 1;
+ } else nago = neighbor->ago;
// set energy/virial flags
// invoke allocate_peratom() if needed for first time
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
if (evflag_atom && !peratom_allocate_flag) {
allocate_peratom();
peratom_allocate_flag = 1;
}
bool success = true;
int flag=PPPM_GPU_API(spread)(nago, atom->nlocal, atom->nlocal +
atom->nghost, atom->x, atom->type, success,
atom->q, domain->boxlo, delxinv, delyinv,
delzinv);
if (!success)
error->one(FLERR,"Insufficient memory on accelerator");
if (flag != 0)
error->one(FLERR,"Out of range atoms - cannot compute PPPM");
// If need per-atom energies/virials, also do particle map on host
// concurrently with GPU calculations
if (evflag_atom) {
memory->destroy(part2grid);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"pppm:part2grid");
particle_map();
}
int i,j;
// convert atoms from box to lamda coords
if (triclinic == 0) boxlo = domain->boxlo;
else {
boxlo = domain->boxlo_lamda;
domain->x2lamda(atom->nlocal);
}
- double t3=MPI_Wtime();
+ double t3 = MPI_Wtime();
// all procs communicate density values from their ghost cells
// to fully sum contribution in their 3d bricks
// remap from 3d decomposition to FFT decomposition
+ cg->reverse_comm(this,REVERSE_RHO);
brick2fft();
// compute potential gradient on my FFT grid and
// portion of e_long on this proc's FFT grid
// return gradients (electric fields) in 3d brick decomposition
poisson();
- // all procs communicate E-field values to fill ghost cells
- // surrounding their 3d bricks
+ // all procs communicate E-field values
+ // to fill ghost cells surrounding their 3d bricks
+
+ if (differentiation_flag == 1) cg->forward_comm(this,FORWARD_AD);
+ else cg->forward_comm(this,FORWARD_IK);
- fillbrick();
+ // extra per-atom energy/virial communication
+
+ if (evflag_atom) {
+ if (differentiation_flag == 1 && vflag_atom)
+ cg_peratom->forward_comm(this,FORWARD_AD_PERATOM);
+ else if (differentiation_flag == 0)
+ cg_peratom->forward_comm(this,FORWARD_IK_PERATOM);
+ }
- poisson_time+=MPI_Wtime()-t3;
+ poisson_time += MPI_Wtime()-t3;
// calculate the force on my particles
FFT_SCALAR qscale = force->qqrd2e * scale;
PPPM_GPU_API(interp)(qscale);
- // Compute per-atom energy/virial on host if requested
+ // per-atom energy/virial
+ // energy includes self-energy correction
if (evflag_atom) {
- fillbrick_peratom();
- fieldforce_peratom();
double *q = atom->q;
int nlocal = atom->nlocal;
if (eflag_atom) {
for (i = 0; i < nlocal; i++) {
eatom[i] *= 0.5;
eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum /
(g_ewald*g_ewald*volume);
eatom[i] *= qscale;
}
}
if (vflag_atom) {
for (i = 0; i < nlocal; i++)
for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale;
}
}
// sum energy across procs and add in volume-dependent term
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qscale;
}
// sum virial across procs
if (vflag_global) {
double virial_all[6];
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
}
// 2d slab correction
if (slabflag) slabcorr();
// convert atoms back from lamda to box coords
if (triclinic) domain->lamda2x(atom->nlocal);
if (kspace_split) PPPM_GPU_API(forces)(atom->f);
}
/* ----------------------------------------------------------------------
- allocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMGPU::allocate()
-{
- memory->create(density_fft,nfft_both,"pppm:density_fft");
- memory->create(greensfn,nfft_both,"pppm:greensfn");
- memory->create(work1,2*nfft_both,"pppm:work1");
- memory->create(work2,2*nfft_both,"pppm:work2");
- memory->create(vg,nfft_both,6,"pppm:vg");
-
- memory->create1d_offset(fkx,nxlo_fft,nxhi_fft,"pppm:fkx");
- memory->create1d_offset(fky,nylo_fft,nyhi_fft,"pppm:fky");
- memory->create1d_offset(fkz,nzlo_fft,nzhi_fft,"pppm:fkz");
-
- memory->create(buf1,nbuf,"pppm:buf1");
- memory->create(buf2,nbuf,"pppm:buf2");
-
- // summation coeffs
-
- memory->create(gf_b,order,"pppm:gf_b");
- memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm:rho1d");
- memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm:drho1d");
- memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm:rho_coeff");
- memory->create2d_offset(drho_coeff,order,(1-order)/2,order/2,
- "pppm:drho_coeff");
-
- // create 2 FFTs and a Remap
- // 1st FFT keeps data in FFT decompostion
- // 2nd FFT returns data in 3d brick decomposition
- // remap takes data from 3d brick to FFT decomposition
-
- int tmp;
-
- fft1 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
- nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
- nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
- 0,0,&tmp);
-
- fft2 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
- nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
- nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
- 0,0,&tmp);
-
- remap = new Remap(lmp,world,
- nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
- nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
- 1,0,0,FFT_PRECISION);
-}
-
-/* ----------------------------------------------------------------------
- deallocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMGPU::deallocate()
-{
- destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
- destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
-
- memory->destroy(density_fft);
- memory->destroy(greensfn);
- memory->destroy(work1);
- memory->destroy(work2);
- memory->destroy(vg);
-
- memory->destroy1d_offset(fkx,nxlo_fft);
- memory->destroy1d_offset(fky,nylo_fft);
- memory->destroy1d_offset(fkz,nzlo_fft);
-
- memory->destroy(buf1);
- memory->destroy(buf2);
-
- memory->destroy(gf_b);
- memory->destroy2d_offset(rho1d,-order/2);
- memory->destroy2d_offset(drho1d,-order/2);
- memory->destroy2d_offset(rho_coeff,(1-order)/2);
- memory->destroy2d_offset(drho_coeff,(1-order)/2);
-
- delete fft1;
- delete fft2;
- delete remap;
-}
-
-
-/* ----------------------------------------------------------------------
- ghost-swap to accumulate full density in brick decomposition
remap density from 3d brick decomposition to FFT decomposition
------------------------------------------------------------------------- */
void PPPMGPU::brick2fft()
{
- int i,n,ix,iy,iz;
- MPI_Request request;
- MPI_Status status;
-
- // pack my ghosts for +x processor
- // pass data to self or +x processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = nxhi_in+1; ix <= nxhi_out; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[0][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = nxlo_in; ix < nxlo_in+nxlo_ghost; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // pack my ghosts for -x processor
- // pass data to self or -x processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = nxlo_out; ix < nxlo_in; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[0][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = nxhi_in-nxhi_ghost+1; ix <= nxhi_in; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // pack my ghosts for +y processor
- // pass data to self or +y processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nyhi_in+1; iy <= nyhi_out; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[1][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // pack my ghosts for -y processor
- // pass data to self or -y processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy < nylo_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[1][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world);
- MPI_Wait(&request,&status);
- }
+ int n,ix,iy,iz;
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // pack my ghosts for +z processor
- // pass data to self or +z processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzhi_in+1; iz <= nzhi_out; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[2][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // pack my ghosts for -z processor
- // pass data to self or -z processor
- // unpack and sum recv data into my real cells
-
- n = 0;
- for (iz = nzlo_out; iz < nzlo_in; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- buf1[n++] = density_brick_gpu[iz][iy][ix];
-
- if (comm->procneigh[2][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = nxlo_in; ix <= nxhi_in; ix++)
- density_brick_gpu[iz][iy][ix] += buf2[n++];
-
- // remap from 3d brick decomposition to FFT decomposition
// copy grabs inner portion of density from 3d brick
// remap could be done as pre-stage of FFT,
// but this works optimally on only double values, not complex values
n = 0;
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
density_fft[n++] = density_brick_gpu[iz][iy][ix];
remap->perform(density_fft,density_fft,work1);
}
-/* ----------------------------------------------------------------------
- Same as base class - needed to call GPU version of fillbrick_.
-------------------------------------------------------------------------- */
-
-void PPPMGPU::fillbrick()
-{
- if (differentiation_flag == 1) fillbrick_ad();
- else fillbrick_ik();
-}
-
-/* ----------------------------------------------------------------------
- ghost-swap to fill ghost cells of my brick with field values
-------------------------------------------------------------------------- */
-
-void PPPMGPU::fillbrick_ik()
-{
- int i,n,ix,iy,iz;
- MPI_Request request;
- MPI_Status status;
-
- // pack my real cells for +z processor
- // pass data to self or +z processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- int x_lo = nxlo_in * 4;
- int x_hi = nxhi_in * 4 + 1;
- for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[2][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz < nzlo_in; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-
- // pack my real cells for -z processor
- // pass data to self or -z processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[2][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzhi_in+1; iz <= nzhi_out; iz++)
- for (iy = nylo_in; iy <= nyhi_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-
- // pack my real cells for +y processor
- // pass data to self or +y processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[1][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy < nylo_in; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-
- // pack my real cells for -y processor
- // pass data to self or -y processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[1][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nyhi_in+1; iy <= nyhi_out; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-
- // pack my real cells for +x processor
- // pass data to self or +x processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- x_lo = (nxhi_in-nxhi_ghost+1) * 4;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[0][1] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- x_lo = nxlo_out * 4;
- x_hi = nxlo_in * 4;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-
- // pack my real cells for -x processor
- // pass data to self or -x processor
- // unpack and sum recv data into my ghost cells
-
- n = 0;
- x_lo = x_hi;
- x_hi = (nxlo_in+nxlo_ghost) * 4;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- buf1[n++] = vd_brick[iz][iy][ix];
- buf1[n++] = vd_brick[iz][iy][ix+1];
- buf1[n++] = vd_brick[iz][iy][ix+2];
- }
-
- if (comm->procneigh[0][0] == me)
- for (i = 0; i < n; i++) buf2[i] = buf1[i];
- else {
- MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world,&request);
- MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world);
- MPI_Wait(&request,&status);
- }
-
- n = 0;
- x_lo = (nxhi_in + 1) * 4;
- x_hi = nxhi_out * 4 + 1;
- for (iz = nzlo_out; iz <= nzhi_out; iz++)
- for (iy = nylo_out; iy <= nyhi_out; iy++)
- for (ix = x_lo; ix < x_hi; ix+=4) {
- vd_brick[iz][iy][ix] = buf2[n++];
- vd_brick[iz][iy][ix+1] = buf2[n++];
- vd_brick[iz][iy][ix+2] = buf2[n++];
- }
-}
-
/* ----------------------------------------------------------------------
Same code as base class - necessary to call GPU version of poisson_ik
------------------------------------------------------------------------- */
void PPPMGPU::poisson()
{
if (differentiation_flag == 1) poisson_ad();
else poisson_ik();
}
/* ----------------------------------------------------------------------
FFT-based Poisson solver
------------------------------------------------------------------------- */
void PPPMGPU::poisson_ik()
{
int i,j,k,n;
double eng;
// transform charge density (r -> k)
n = 0;
for (i = 0; i < nfft; i++) {
work1[n++] = density_fft[i];
work1[n++] = ZEROF;
}
fft1->compute(work1,work1,1);
// if requested, compute energy and virial contribution
double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
double s2 = scaleinv*scaleinv;
if (eflag_global || vflag_global) {
if (vflag_global) {
n = 0;
for (i = 0; i < nfft; i++) {
eng = s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]);
for (j = 0; j < 6; j++) virial[j] += eng*vg[i][j];
if (eflag_global) energy += eng;
n += 2;
}
} else {
n = 0;
for (i = 0; i < nfft; i++) {
energy +=
s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]);
n += 2;
}
}
}
// scale by 1/total-grid-pts to get rho(k)
// multiply by Green's function to get V(k)
n = 0;
for (i = 0; i < nfft; i++) {
work1[n++] *= scaleinv * greensfn[i];
work1[n++] *= scaleinv * greensfn[i];
}
// extra FFTs for per-atom energy/virial
if (evflag_atom) poisson_peratom();
// compute gradients of V(r) in each of 3 dims by transformimg -ik*V(k)
// FFT leaves data in 3d brick decomposition
// copy it into inner portion of vdx,vdy,vdz arrays
// x direction gradient
n = 0;
for (k = nzlo_fft; k <= nzhi_fft; k++)
for (j = nylo_fft; j <= nyhi_fft; j++)
for (i = nxlo_fft; i <= nxhi_fft; i++) {
work2[n] = fkx[i]*work1[n+1];
work2[n+1] = -fkx[i]*work1[n];
n += 2;
}
fft2->compute(work2,work2,-1);
n = 0;
int x_hi = nxhi_in * 4 + 3;
for (k = nzlo_in; k <= nzhi_in; k++)
for (j = nylo_in; j <= nyhi_in; j++)
for (i = nxlo_in * 4; i < x_hi; i+=4) {
vd_brick[k][j][i] = work2[n];
n += 2;
}
// y direction gradient
n = 0;
for (k = nzlo_fft; k <= nzhi_fft; k++)
for (j = nylo_fft; j <= nyhi_fft; j++)
for (i = nxlo_fft; i <= nxhi_fft; i++) {
work2[n] = fky[j]*work1[n+1];
work2[n+1] = -fky[j]*work1[n];
n += 2;
}
fft2->compute(work2,work2,-1);
n = 0;
for (k = nzlo_in; k <= nzhi_in; k++)
for (j = nylo_in; j <= nyhi_in; j++)
for (i = nxlo_in * 4 + 1; i < x_hi; i+=4) {
vd_brick[k][j][i] = work2[n];
n += 2;
}
// z direction gradient
n = 0;
for (k = nzlo_fft; k <= nzhi_fft; k++)
for (j = nylo_fft; j <= nyhi_fft; j++)
for (i = nxlo_fft; i <= nxhi_fft; i++) {
work2[n] = fkz[k]*work1[n+1];
work2[n+1] = -fkz[k]*work1[n];
n += 2;
}
fft2->compute(work2,work2,-1);
n = 0;
for (k = nzlo_in; k <= nzhi_in; k++)
for (j = nylo_in; j <= nyhi_in; j++)
for (i = nxlo_in * 4 + 2; i < x_hi; i+=4) {
vd_brick[k][j][i] = work2[n];
n += 2;
}
}
/* ----------------------------------------------------------------------
- Create array using offsets from pinned memory allocation
+ pack own values to buf to send to another proc
+------------------------------------------------------------------------- */
+
+void PPPMGPU::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+ int n = 0;
+
+ if (flag == FORWARD_IK) {
+ int offset;
+ FFT_SCALAR *src = &vd_brick[nzlo_out][nylo_out][4*nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ offset = 4*list[i];
+ buf[n++] = src[offset++];
+ buf[n++] = src[offset++];
+ buf[n++] = src[offset];
+ }
+ } else if (flag == FORWARD_AD) {
+ FFT_SCALAR *src = &u_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++)
+ buf[i] = src[list[i]];
+ } else if (flag == FORWARD_IK_PERATOM) {
+ FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ if (eflag_atom) buf[n++] = esrc[list[i]];
+ if (vflag_atom) {
+ buf[n++] = v0src[list[i]];
+ buf[n++] = v1src[list[i]];
+ buf[n++] = v2src[list[i]];
+ buf[n++] = v3src[list[i]];
+ buf[n++] = v4src[list[i]];
+ buf[n++] = v5src[list[i]];
+ }
+ }
+ } else if (flag == FORWARD_AD_PERATOM) {
+ FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ buf[n++] = v0src[list[i]];
+ buf[n++] = v1src[list[i]];
+ buf[n++] = v2src[list[i]];
+ buf[n++] = v3src[list[i]];
+ buf[n++] = v4src[list[i]];
+ buf[n++] = v5src[list[i]];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ unpack another proc's own values from buf and set own ghost values
+------------------------------------------------------------------------- */
+
+void PPPMGPU::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+ int n = 0;
+
+ if (flag == FORWARD_IK) {
+ int offset;
+ FFT_SCALAR *dest = &vdx_brick[nzlo_out][nylo_out][4*nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ offset = 4*list[i];
+ dest[offset++] = buf[n++];
+ dest[offset++] = buf[n++];
+ dest[offset] = buf[n++];
+ }
+ } else if (flag == FORWARD_AD) {
+ FFT_SCALAR *dest = &u_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++)
+ dest[list[i]] = buf[n++];
+ } else if (flag == FORWARD_IK_PERATOM) {
+ FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ if (eflag_atom) esrc[list[i]] = buf[n++];
+ if (vflag_atom) {
+ v0src[list[i]] = buf[n++];
+ v1src[list[i]] = buf[n++];
+ v2src[list[i]] = buf[n++];
+ v3src[list[i]] = buf[n++];
+ v4src[list[i]] = buf[n++];
+ v5src[list[i]] = buf[n++];
+ }
+ }
+ } else if (flag == FORWARD_AD_PERATOM) {
+ FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+ FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++) {
+ v0src[list[i]] = buf[n++];
+ v1src[list[i]] = buf[n++];
+ v2src[list[i]] = buf[n++];
+ v3src[list[i]] = buf[n++];
+ v4src[list[i]] = buf[n++];
+ v5src[list[i]] = buf[n++];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ pack ghost values into buf to send to another proc
+------------------------------------------------------------------------- */
+
+void PPPMGPU::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+ if (flag == REVERSE_RHO) {
+ FFT_SCALAR *src = &density_brick_gpu[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++)
+ buf[i] = src[list[i]];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ unpack another proc's ghost values from buf and add to own values
+------------------------------------------------------------------------- */
+
+void PPPMGPU::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+ if (flag == REVERSE_RHO) {
+ FFT_SCALAR *dest = &density_brick_gpu[nzlo_out][nylo_out][nxlo_out];
+ for (int i = 0; i < nlist; i++)
+ dest[list[i]] += buf[i];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ create array using offsets from pinned memory allocation
------------------------------------------------------------------------- */
FFT_SCALAR ***PPPMGPU::create_3d_offset(int n1lo, int n1hi, int n2lo, int n2hi,
- int n3lo, int n3hi, const char *name,
- FFT_SCALAR *data, int vec_length)
+ int n3lo, int n3hi, const char *name,
+ FFT_SCALAR *data, int vec_length)
{
int i,j;
int n1 = n1hi - n1lo + 1;
int n2 = n2hi - n2lo + 1;
int n3 = n3hi - n3lo + 1;
FFT_SCALAR **plane = (FFT_SCALAR **)
memory->smalloc(n1*n2*sizeof(FFT_SCALAR *),name);
FFT_SCALAR ***array = (FFT_SCALAR ***)
memory->smalloc(n1*sizeof(FFT_SCALAR **),name);
int n = 0;
for (i = 0; i < n1; i++) {
array[i] = &plane[i*n2];
for (j = 0; j < n2; j++) {
plane[i*n2+j] = &data[n];
n += n3*vec_length;
}
}
for (i = 0; i < n1*n2; i++) array[0][i] -= n3lo*vec_length;
for (i = 0; i < n1; i++) array[i] -= n2lo;
return array-n1lo;
}
/* ----------------------------------------------------------------------
3d memory offsets
------------------------------------------------------------------------- */
void PPPMGPU::destroy_3d_offset(FFT_SCALAR ***array, int n1_offset,
int n2_offset)
{
if (array == NULL) return;
memory->sfree(&array[n1_offset][n2_offset]);
memory->sfree(array + n1_offset);
}
/* ----------------------------------------------------------------------
memory usage of local arrays
------------------------------------------------------------------------- */
double PPPMGPU::memory_usage()
{
- double bytes = nmax*3 * sizeof(double);
- int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
- (nzhi_out-nzlo_out+1);
- bytes += 4 * nbrick * sizeof(FFT_SCALAR);
- bytes += 6 * nfft_both * sizeof(double);
- bytes += nfft_both * sizeof(double);
- bytes += nfft_both*5 * sizeof(FFT_SCALAR);
- bytes += 2 * nbuf * sizeof(double);
-
- if (peratom_allocate_flag) {
- bytes += 7 * nbrick * sizeof(FFT_SCALAR);
- bytes += 2 * nbuf_peratom * sizeof(FFT_SCALAR);
- }
+ double bytes = PPPM::memory_usage();
+
+ // NOTE: add tallying here for density_brick_gpu and vd_brick
+ // could subtract out density_brick and vdxyz_brick if freed them above
+ // it the net efffect is zero, do nothing
return bytes + PPPM_GPU_API(bytes)();
}
/* ----------------------------------------------------------------------
perform and time the 4 FFTs required for N timesteps
------------------------------------------------------------------------- */
int PPPMGPU::timing(int n, double &time3d, double &time1d) {
if (im_real_space) {
time3d = 1.0;
time1d = 1.0;
return 4;
}
PPPM::timing(n,time3d,time1d);
return 4;
}
/* ----------------------------------------------------------------------
adjust PPPM coeffs, called initially and whenever volume has changed
------------------------------------------------------------------------- */
void PPPMGPU::setup()
{
if (im_real_space) return;
PPPM::setup();
}
diff --git a/src/GPU/pppm_gpu.h b/src/GPU/pppm_gpu.h
index 4c9e8441a..98d62ffa2 100644
--- a/src/GPU/pppm_gpu.h
+++ b/src/GPU/pppm_gpu.h
@@ -1,100 +1,99 @@
/* ----------------------------------------------------------------------
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 KSPACE_CLASS
KSpaceStyle(pppm/gpu,PPPMGPU)
#else
#ifndef LMP_PPPM_GPU_H
#define LMP_PPPM_GPU_H
#include "pppm.h"
namespace LAMMPS_NS {
class PPPMGPU : public PPPM {
public:
PPPMGPU(class LAMMPS *, int, char **);
virtual ~PPPMGPU();
void init();
void setup();
void compute(int, int);
int timing(int, double &, double &);
double memory_usage();
protected:
-
FFT_SCALAR ***density_brick_gpu, ***vd_brick;
bool kspace_split, im_real_space;
+ int old_nlocal;
+ double poisson_time;
- void allocate();
- void deallocate();
void brick2fft();
- void fillbrick();
- void fillbrick_ik();
void poisson();
void poisson_ik();
- int old_nlocal;
- double poisson_time;
+ void pack_forward(int, FFT_SCALAR *, int, int *);
+ void unpack_forward(int, FFT_SCALAR *, int, int *);
+ void pack_reverse(int, FFT_SCALAR *, int, int *);
+ void unpack_reverse(int, FFT_SCALAR *, int, int *);
FFT_SCALAR ***create_3d_offset(int, int, int, int, int, int, const char *,
FFT_SCALAR *, int);
void destroy_3d_offset(FFT_SCALAR ***, int, int);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot (yet) do analytic differentiation with pppm/gpu.
Self-explanatory.
E: Cannot use order greater than 8 with pppm/gpu.
Self-explanatory.
E: Insufficient memory on accelerator
There is insufficient memory on one of the devices specified for the gpu
package
E: Out of range atoms - cannot compute PPPM
One or more atoms are attempting to map their charge to a PPPM grid
point that is not owned by a processor. This is likely for one of two
reasons, both of them bad. First, it may mean that an atom near the
boundary of a processor's sub-domain has moved more than 1/2 the
"neighbor skin distance"_neighbor.html without neighbor lists being
rebuilt and atoms being migrated to new processors. This also means
you may be missing pairwise interactions that need to be computed.
The solution is to change the re-neighboring criteria via the
"neigh_modify"_neigh_modify command. The safest settings are "delay 0
every 1 check yes". Second, it may mean that an atom has moved far
outside a processor's sub-domain or even the entire simulation box.
This indicates bad physics, e.g. due to highly overlapping atoms, too
large a timestep, etc.
*/
diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h
index 2f5dc2241..a7d0dfb2c 100644
--- a/src/KSPACE/pppm.h
+++ b/src/KSPACE/pppm.h
@@ -1,306 +1,306 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifdef KSPACE_CLASS
KSpaceStyle(pppm,PPPM)
#else
#ifndef LMP_PPPM_H
#define LMP_PPPM_H
#include "lmptype.h"
#include "mpi.h"
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
#include "kspace.h"
namespace LAMMPS_NS {
class PPPM : public KSpace {
public:
PPPM(class LAMMPS *, int, char **);
virtual ~PPPM();
virtual void init();
virtual void setup();
void setup_grid();
virtual void compute(int, int);
virtual int timing_1d(int, double &);
virtual int timing_3d(int, double &);
virtual double memory_usage();
virtual void compute_group_group(int, int, int);
protected:
int me,nprocs;
int nfactors;
int *factors;
double qsum,qsqsum,q2;
double cutoff;
double volume;
double delxinv,delyinv,delzinv,delvolinv;
double shift,shiftone;
int peratom_allocate_flag;
int nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in;
int nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out;
int nxlo_ghost,nxhi_ghost,nylo_ghost,nyhi_ghost,nzlo_ghost,nzhi_ghost;
int nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft;
int nlower,nupper;
int ngrid,nfft,nfft_both;
FFT_SCALAR ***density_brick;
FFT_SCALAR ***vdx_brick,***vdy_brick,***vdz_brick;
FFT_SCALAR ***u_brick;
FFT_SCALAR ***v0_brick,***v1_brick,***v2_brick;
FFT_SCALAR ***v3_brick,***v4_brick,***v5_brick;
double *greensfn;
double **vg;
double *fkx,*fky,*fkz;
FFT_SCALAR *density_fft;
FFT_SCALAR *work1,*work2;
double *gf_b;
FFT_SCALAR **rho1d,**rho_coeff,**drho1d,**drho_coeff;
double *sf_precoeff1, *sf_precoeff2, *sf_precoeff3;
double *sf_precoeff4, *sf_precoeff5, *sf_precoeff6;
double sf_coeff[6]; // coefficients for calculating ad self-forces
double **acons;
// group-group interactions
int group_allocate_flag;
FFT_SCALAR ***density_A_brick,***density_B_brick;
FFT_SCALAR *density_A_fft,*density_B_fft;
class FFT3d *fft1,*fft2;
class Remap *remap;
class CommGrid *cg;
class CommGrid *cg_peratom;
int **part2grid; // storage for particle -> grid mapping
int nmax;
int triclinic; // domain settings, orthog or triclinic
double *boxlo;
// TIP4P settings
int typeH,typeO; // atom types of TIP4P water H and O atoms
double qdist; // distance from O site to negative charge
double alpha; // geometric factor
void set_grid_global();
void set_grid_local();
void adjust_gewald();
double newton_raphson_f();
double derivf();
double final_accuracy();
virtual void allocate();
virtual void allocate_peratom();
virtual void deallocate();
virtual void deallocate_peratom();
int factorable(int);
double compute_df_kspace();
double estimate_ik_error(double, double, bigint);
double compute_qopt();
void compute_gf_denom();
void compute_gf_ik();
void compute_gf_ad();
void compute_sf_precoeff();
virtual void particle_map();
virtual void make_rho();
virtual void brick2fft();
virtual void poisson();
void poisson_ik();
void poisson_ad();
virtual void fieldforce();
void fieldforce_ik();
void fieldforce_ad();
virtual void poisson_peratom();
virtual void fieldforce_peratom();
void procs2grid2d(int,int,int,int *, int*);
void compute_rho1d(const FFT_SCALAR &, const FFT_SCALAR &,
const FFT_SCALAR &);
void compute_drho1d(const FFT_SCALAR &, const FFT_SCALAR &,
const FFT_SCALAR &);
void compute_rho_coeff();
void slabcorr();
// grid communication
- void pack_forward(int, FFT_SCALAR *, int, int *);
- void unpack_forward(int, FFT_SCALAR *, int, int *);
- void pack_reverse(int, FFT_SCALAR *, int, int *);
- void unpack_reverse(int, FFT_SCALAR *, int, int *);
+ virtual void pack_forward(int, FFT_SCALAR *, int, int *);
+ virtual void unpack_forward(int, FFT_SCALAR *, int, int *);
+ virtual void pack_reverse(int, FFT_SCALAR *, int, int *);
+ virtual void unpack_reverse(int, FFT_SCALAR *, int, int *);
// group-group interactions
virtual void allocate_groups();
virtual void deallocate_groups();
virtual void make_rho_groups(int, int, int);
virtual void poisson_groups(int);
/* ----------------------------------------------------------------------
denominator for Hockney-Eastwood Green's function
of x,y,z = sin(kx*deltax/2), etc
inf n-1
S(n,k) = Sum W(k+pi*j)**2 = Sum b(l)*(z*z)**l
j=-inf l=0
= -(z*z)**n /(2n-1)! * (d/dx)**(2n-1) cot(x) at z = sin(x)
gf_b = denominator expansion coeffs
------------------------------------------------------------------------- */
inline double gf_denom(const double &x, const double &y,
const double &z) const {
double sx,sy,sz;
sz = sy = sx = 0.0;
for (int l = order-1; l >= 0; l--) {
sx = gf_b[l] + sx*x;
sy = gf_b[l] + sy*y;
sz = gf_b[l] + sz*z;
}
double s = sx*sy*sz;
return s*s;
};
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot (yet) use PPPM with triclinic box
This feature is not yet supported.
E: Cannot use PPPM with 2d simulation
The kspace style pppm cannot be used in 2d simulations. You can use
2d PPPM in a 3d simulation; see the kspace_modify command.
E: Kspace style requires atom attribute q
The atom style defined does not have these attributes.
E: Cannot use nonperiodic boundaries with PPPM
For kspace style pppm, all 3 dimensions must have periodic boundaries
unless you use the kspace_modify command to define a 2d slab with a
non-periodic z dimension.
E: Incorrect boundaries with slab PPPM
Must have periodic x,y dimensions and non-periodic z dimension to use
2d slab option with PPPM.
E: PPPM order cannot be < 2 or > than %d
This is a limitation of the PPPM implementation in LAMMPS.
E: KSpace style is incompatible with Pair style
Setting a kspace style requires that a pair style with a long-range
Coulombic component be selected.
E: Bond and angle potentials must be defined for TIP4P
Cannot use TIP4P pair potential unless bond and angle potentials
are defined.
E: Bad TIP4P angle type for PPPM/TIP4P
Specified angle type is not valid.
E: Bad TIP4P bond type for PPPM/TIP4P
Specified bond type is not valid.
E: Cannot use kspace solver on system with no charge
No atoms in system have a non-zero charge.
W: System is not charge neutral, net charge = %g
The total charge on all atoms on the system is not 0.0, which
is not valid for Ewald or PPPM.
W: Reducing PPPM order b/c stencil extends beyond neighbor processor
LAMMPS is attempting this in order to allow the simulation
to run. It should not effect the PPPM accuracy.
E: PPPM grid is too large
The global PPPM grid is larger than OFFSET in one or more dimensions.
OFFSET is currently set to 4096. You likely need to decrease the
requested accuracy.
E: PPPM order has been reduced to 0
LAMMPS has attempted to reduce the PPPM order to enable the simulation
to run, but can reduce the order no further. Try increasing the
accuracy of PPPM by reducing the tolerance size, thus inducing a
larger PPPM grid.
E: KSpace accuracy must be > 0
The kspace accuracy designated in the input must be greater than zero.
E: Cannot compute PPPM G
LAMMPS failed to compute a valid approximation for the PPPM g_ewald
factor that partitions the computation between real space and k-space.
E: Out of range atoms - cannot compute PPPM
One or more atoms are attempting to map their charge to a PPPM grid
point that is not owned by a processor. This is likely for one of two
reasons, both of them bad. First, it may mean that an atom near the
boundary of a processor's sub-domain has moved more than 1/2 the
"neighbor skin distance"_neighbor.html without neighbor lists being
rebuilt and atoms being migrated to new processors. This also means
you may be missing pairwise interactions that need to be computed.
The solution is to change the re-neighboring criteria via the
"neigh_modify"_neigh_modify command. The safest settings are "delay 0
every 1 check yes". Second, it may mean that an atom has moved far
outside a processor's sub-domain or even the entire simulation box.
This indicates bad physics, e.g. due to highly overlapping atoms, too
large a timestep, etc.
E: Cannot (yet) use K-space slab correction with compute group/group
This option is not yet supported.
*/
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 658ac11ca..d14d5084a 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -1,1576 +1,1576 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "ctype.h"
#include "read_data.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "error.h"
#include "memory.h"
#include "special.h"
using namespace LAMMPS_NS;
#define MAXLINE 256
#define LB_FACTOR 1.1
#define CHUNK 1024
#define DELTA 4 // must be 2 or larger
// customize for new sections
#define NSECTIONS 23 // change when add to header::section_keywords
/* ---------------------------------------------------------------------- */
ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
line = new char[MAXLINE];
keyword = new char[MAXLINE];
buffer = new char[CHUNK*MAXLINE];
narg = maxarg = 0;
arg = NULL;
// customize for new sections
// pointers to atom styles that store extra info
nellipsoids = 0;
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
nlines = 0;
avec_line = (AtomVecLine *) atom->style_match("line");
ntris = 0;
avec_tri = (AtomVecTri *) atom->style_match("tri");
}
/* ---------------------------------------------------------------------- */
ReadData::~ReadData()
{
delete [] line;
delete [] keyword;
delete [] buffer;
memory->sfree(arg);
}
/* ---------------------------------------------------------------------- */
void ReadData::command(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Illegal read_data command");
if (domain->box_exist)
error->all(FLERR,"Cannot read_data after simulation box is defined");
if (domain->dimension == 2 && domain->zperiodic == 0)
error->all(FLERR,"Cannot run 2d simulation with nonperiodic Z dimension");
// fixes that process data file info
nfix = 0;
fix_index = NULL;
fix_header = NULL;
fix_section = NULL;
int iarg = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"fix") == 0) {
if (iarg+4 > narg)
error->all(FLERR,"Illegal read_data command");
memory->grow(fix_index,nfix+1,"read_data:fix_index");
fix_header = (char **)
memory->srealloc(fix_header,(nfix+1)*sizeof(char *),
"read_data:fix_header");
fix_section = (char **)
memory->srealloc(fix_section,(nfix+1)*sizeof(char *),
"read_data:fix_section");
fix_index[nfix] = modify->find_fix(arg[iarg+1]);
if (fix_index[nfix] < 0)
error->all(FLERR,"Fix ID for read_data does not exist");
int n = strlen(arg[iarg+2]) + 1;
fix_header[nfix] = new char[n];
strcpy(fix_header[nfix],arg[iarg+2]);
n = strlen(arg[iarg+3]) + 1;
fix_section[nfix] = new char[n];
strcpy(fix_section[nfix],arg[iarg+3]);
nfix++;
iarg += 4;
} else error->all(FLERR,"Illegal read_data command");
}
// scan data file to determine max topology needed per atom
// allocate initial topology arrays
if (atom->molecular) {
if (me == 0) {
if (screen) fprintf(screen,"Scanning data file ...\n");
open(arg[0]);
header(0);
scan(atom->bond_per_atom,atom->angle_per_atom,
atom->dihedral_per_atom,atom->improper_per_atom);
if (compressed) pclose(fp);
else fclose(fp);
atom->bond_per_atom += atom->extra_bond_per_atom;
}
MPI_Bcast(&atom->bond_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->angle_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->dihedral_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->improper_per_atom,1,MPI_INT,0,world);
} else
atom->bond_per_atom = atom->angle_per_atom =
atom->dihedral_per_atom = atom->improper_per_atom = 0;
// read header info
if (me == 0) {
if (screen) fprintf(screen,"Reading data file ...\n");
open(arg[0]);
}
header(1);
domain->box_exist = 1;
// problem setup using info from header
update->ntimestep = 0;
int n;
if (comm->nprocs == 1) n = static_cast<int> (atom->natoms);
else n = static_cast<int> (LB_FACTOR * atom->natoms / comm->nprocs);
atom->allocate_type_arrays();
atom->avec->grow(n);
n = atom->nmax;
domain->print_box(" ");
domain->set_initial_box();
domain->set_global_box();
comm->set_proc_grid();
domain->set_local_box();
// customize for new sections
// read rest of file in free format
int atomflag = 0;
while (strlen(keyword)) {
// allow special fixes first chance to match and process the section
// if fix matches, continue to next section
if (nfix) {
for (n = 0; n < nfix; n++)
if (strstr(line,fix_section[n])) {
bigint nlines =
modify->fix[fix_index[n]]->read_data_skip_lines(keyword);
fix(n,keyword,nlines);
parse_keyword(0,1);
break;
}
if (n < nfix) continue;
}
if (strcmp(keyword,"Atoms") == 0) {
atoms();
atomflag = 1;
} else if (strcmp(keyword,"Velocities") == 0) {
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Velocities");
velocities();
} else if (strcmp(keyword,"Ellipsoids") == 0) {
if (!avec_ellipsoid)
error->all(FLERR,"Invalid data file section: Ellipsoids");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Ellipsoids");
bonus(nellipsoids,(AtomVec *) avec_ellipsoid,"ellipsoids");
} else if (strcmp(keyword,"Lines") == 0) {
if (!avec_line)
error->all(FLERR,"Invalid data file section: Lines");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Lines");
bonus(nlines,(AtomVec *) avec_line,"lines");
} else if (strcmp(keyword,"Triangles") == 0) {
if (!avec_tri)
error->all(FLERR,"Invalid data file section: Triangles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles");
bonus(ntris,(AtomVec *) avec_tri,"triangles");
} else if (strcmp(keyword,"Bonds") == 0) {
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Invalid data file section: Bonds");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bonds");
bonds();
} else if (strcmp(keyword,"Angles") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: Angles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Angles");
angles();
} else if (strcmp(keyword,"Dihedrals") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: Dihedrals");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Dihedrals");
dihedrals();
} else if (strcmp(keyword,"Impropers") == 0) {
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Invalid data file section: Impropers");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Impropers");
impropers();
} else if (strcmp(keyword,"Masses") == 0) {
mass();
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
error->all(FLERR,"Must define pair_style before Pair Coeffs");
paircoeffs();
} else if (strcmp(keyword,"Bond Coeffs") == 0) {
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Invalid data file section: Bond Coeffs");
if (force->bond == NULL)
error->all(FLERR,"Must define bond_style before Bond Coeffs");
bondcoeffs();
} else if (strcmp(keyword,"Angle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: Angle Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before Angle Coeffs");
anglecoeffs(0);
} else if (strcmp(keyword,"Dihedral Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: Dihedral Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before Dihedral Coeffs");
dihedralcoeffs(0);
} else if (strcmp(keyword,"Improper Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Invalid data file section: Improper Coeffs");
if (force->improper == NULL)
error->all(FLERR,"Must define improper_style before Improper Coeffs");
impropercoeffs(0);
} else if (strcmp(keyword,"BondBond Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: BondBond Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before BondBond Coeffs");
anglecoeffs(1);
} else if (strcmp(keyword,"BondAngle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: BondAngle Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before BondAngle Coeffs");
anglecoeffs(2);
} else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
dihedralcoeffs(1);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
dihedralcoeffs(2);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
dihedralcoeffs(3);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
dihedralcoeffs(4);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: BondBond13 Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before BondBond13 Coeffs");
dihedralcoeffs(5);
} else if (strcmp(keyword,"AngleAngle Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngle Coeffs");
if (force->improper == NULL)
error->all(FLERR,"Must define improper_style before AngleAngle Coeffs");
impropercoeffs(1);
} else {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->all(FLERR,str);
}
parse_keyword(0,1);
}
// close file
if (me == 0) {
if (compressed) pclose(fp);
else fclose(fp);
}
// error if natoms > 0 yet no atoms were read
if (atom->natoms > 0 && atomflag == 0)
error->all(FLERR,"No atoms in data file");
// create bond topology now that system is defined
if (atom->molecular) {
Special special(lmp);
special.build();
}
}
/* ----------------------------------------------------------------------
read free-format header of data file
if flag = 0, only called by proc 0
if flag = 1, called by all procs so bcast lines as read them
1st line and blank lines are skipped
non-blank lines are checked for header keywords and leading value is read
header ends with EOF or non-blank line containing no header keyword
if EOF, line is set to blank line
else line has first keyword line for rest of file
------------------------------------------------------------------------- */
void ReadData::header(int flag)
{
int n;
char *ptr;
// customize for new sections
const char *section_keywords[NSECTIONS] =
{"Atoms","Velocities","Ellipsoids","Lines","Triangles",
"Bonds","Angles","Dihedrals","Impropers",
"Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
"BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs",
"EndBondTorsion Coeffs","AngleTorsion Coeffs",
"AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"};
// skip 1st line of file
if (me == 0) {
char *eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
}
// customize for new header lines
while (1) {
// read a line and bcast length if flag is set
if (me == 0) {
if (fgets(line,MAXLINE,fp) == NULL) n = 0;
else n = strlen(line) + 1;
}
if (flag) MPI_Bcast(&n,1,MPI_INT,0,world);
// if n = 0 then end-of-file so return with blank line
if (n == 0) {
line[0] = '\0';
return;
}
// bcast line if flag is set
if (flag) MPI_Bcast(line,n,MPI_CHAR,0,world);
// trim anything from '#' onward
// if line is blank, continue
if (ptr = strchr(line,'#')) *ptr = '\0';
if (strspn(line," \t\n\r") == strlen(line)) continue;
// allow special fixes first chance to match and process the line
// if fix matches, continue to next header line
if (nfix) {
for (n = 0; n < nfix; n++)
if (strstr(line,fix_header[n])) {
modify->fix[fix_index[n]]->read_data_header(line);
break;
}
if (n < nfix) continue;
}
// search line for header keyword and set corresponding variable
if (strstr(line,"atoms")) sscanf(line,BIGINT_FORMAT,&atom->natoms);
// check for these first
// otherwise "triangles" will be matched as "angles"
else if (strstr(line,"ellipsoids")) {
if (!avec_ellipsoid && me == 0)
error->one(FLERR,"No ellipsoids allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nellipsoids);
} else if (strstr(line,"lines")) {
if (!avec_line && me == 0)
error->one(FLERR,"No lines allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nlines);
} else if (strstr(line,"triangles")) {
if (!avec_tri && me == 0)
error->one(FLERR,"No triangles allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&ntris);
}
else if (strstr(line,"bonds")) sscanf(line,BIGINT_FORMAT,&atom->nbonds);
else if (strstr(line,"angles")) sscanf(line,BIGINT_FORMAT,&atom->nangles);
else if (strstr(line,"dihedrals")) sscanf(line,BIGINT_FORMAT,
&atom->ndihedrals);
else if (strstr(line,"impropers")) sscanf(line,BIGINT_FORMAT,
&atom->nimpropers);
else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes);
else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes);
else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes);
else if (strstr(line,"dihedral types"))
sscanf(line,"%d",&atom->ndihedraltypes);
else if (strstr(line,"improper types"))
sscanf(line,"%d",&atom->nimpropertypes);
else if (strstr(line,"extra bond per atom"))
sscanf(line,"%d",&atom->extra_bond_per_atom);
else if (strstr(line,"xlo xhi"))
sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]);
else if (strstr(line,"ylo yhi"))
sscanf(line,"%lg %lg",&domain->boxlo[1],&domain->boxhi[1]);
else if (strstr(line,"zlo zhi"))
sscanf(line,"%lg %lg",&domain->boxlo[2],&domain->boxhi[2]);
else if (strstr(line,"xy xz yz")) {
domain->triclinic = 1;
sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz);
} else break;
}
// error check on total system size
if (atom->natoms < 0 || atom->natoms > MAXBIGINT ||
atom->nbonds < 0 || atom->nbonds > MAXBIGINT ||
atom->nangles < 0 || atom->nangles > MAXBIGINT ||
atom->ndihedrals < 0 || atom->ndihedrals > MAXBIGINT ||
atom->nimpropers < 0 || atom->nimpropers > MAXBIGINT) {
if (me == 0) error->one(FLERR,"System in data file is too big");
}
// check that exiting string is a valid section keyword
parse_keyword(1,flag);
for (n = 0; n < NSECTIONS; n++)
if (strcmp(keyword,section_keywords[n]) == 0) break;
if (n == NSECTIONS && me == 0) {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->one(FLERR,str);
}
// error check on consistency of header values
if ((atom->nbonds || atom->nbondtypes) &&
atom->avec->bonds_allow == 0 && me == 0)
error->one(FLERR,"No bonds allowed with this atom style");
if ((atom->nangles || atom->nangletypes) &&
atom->avec->angles_allow == 0 && me == 0)
error->one(FLERR,"No angles allowed with this atom style");
if ((atom->ndihedrals || atom->ndihedraltypes) &&
atom->avec->dihedrals_allow == 0 && me == 0)
error->one(FLERR,"No dihedrals allowed with this atom style");
if ((atom->nimpropers || atom->nimpropertypes) &&
atom->avec->impropers_allow == 0 && me == 0)
error->one(FLERR,"No impropers allowed with this atom style");
if (atom->nbonds > 0 && atom->nbondtypes <= 0 && me == 0)
error->one(FLERR,"Bonds defined but no bond types");
if (atom->nangles > 0 && atom->nangletypes <= 0 && me == 0)
error->one(FLERR,"Angles defined but no angle types");
if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0 && me == 0)
error->one(FLERR,"Dihedrals defined but no dihedral types");
if (atom->nimpropers > 0 && atom->nimpropertypes <= 0 && me == 0)
error->one(FLERR,"Impropers defined but no improper types");
}
/* ----------------------------------------------------------------------
read all atoms
------------------------------------------------------------------------- */
void ReadData::atoms()
{
int i,m,nchunk;
bigint nread = 0;
bigint natoms = atom->natoms;
while (nread < natoms) {
if (natoms-nread > CHUNK) nchunk = CHUNK;
else nchunk = natoms-nread;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_atoms(nchunk,buffer);
nread += nchunk;
}
// check that all atoms were assigned correctly
bigint tmp = atom->nlocal;
MPI_Allreduce(&tmp,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms);
}
if (natoms != atom->natoms)
error->all(FLERR,"Did not assign all atoms correctly");
// if any atom ID < 0, error
// if all atom IDs = 0, tag_enable = 0
// if any atom ID > 0, error if any atom ID == 0
// not checking if atom IDs > natoms or are unique
int nlocal = atom->nlocal;
int *tag = atom->tag;
int flag = 0;
for (int i = 0; i < nlocal; i++)
if (tag[i] < 0) flag = 1;
int flag_all;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
if (flag_all)
error->all(FLERR,"Invalid atom ID in Atoms section of data file");
flag = 0;
for (int i = 0; i < nlocal; i++)
if (tag[i] > 0) flag = 1;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
if (flag_all == 0) atom->tag_enable = 0;
if (atom->tag_enable) {
flag = 0;
for (int i = 0; i < nlocal; i++)
if (tag[i] == 0) flag = 1;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
if (flag_all)
error->all(FLERR,"Invalid atom ID in Atoms section of data file");
}
// create global mapping
if (atom->map_style) {
atom->map_init();
atom->map_set();
}
}
/* ----------------------------------------------------------------------
read all velocities
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
void ReadData::velocities()
{
int i,m,nchunk;
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_style = 1;
atom->map_init();
atom->map_set();
}
bigint nread = 0;
bigint natoms = atom->natoms;
while (nread < natoms) {
if (natoms-nread > CHUNK) nchunk = CHUNK;
else nchunk = natoms-nread;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_vels(nchunk,buffer);
nread += nchunk;
}
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " velocities\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " velocities\n",natoms);
}
}
/* ----------------------------------------------------------------------
read all bonus data
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
{
int i,m,nchunk;
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_style = 1;
atom->map_init();
atom->map_set();
}
bigint nread = 0;
bigint natoms = nbonus;
while (nread < natoms) {
if (natoms-nread > CHUNK) nchunk = CHUNK;
else nchunk = natoms-nread;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_bonus(nchunk,buffer,ptr);
nread += nchunk;
}
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " %s\n",natoms,type);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " %s\n",natoms,type);
}
}
/* ---------------------------------------------------------------------- */
void ReadData::bonds()
{
int i,m,nchunk;
bigint nread = 0;
bigint nbonds = atom->nbonds;
while (nread < nbonds) {
nchunk = MIN(nbonds-nread,CHUNK);
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_bonds(nchunk,buffer);
nread += nchunk;
}
// check that bonds were assigned correctly
int nlocal = atom->nlocal;
bigint sum;
bigint n = 0;
for (i = 0; i < nlocal; i++) n += atom->num_bond[i];
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 2;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor);
}
if (sum != factor*atom->nbonds)
error->all(FLERR,"Bonds assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
void ReadData::angles()
{
int i,m,nchunk;
bigint nread = 0;
bigint nangles = atom->nangles;
while (nread < nangles) {
nchunk = MIN(nangles-nread,CHUNK);
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_angles(nchunk,buffer);
nread += nchunk;
}
// check that ang
int nlocal = atom->nlocal;
bigint sum;
bigint n = 0;
for (i = 0; i < nlocal; i++) n += atom->num_angle[i];
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 3;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor);
}
if (sum != factor*atom->nangles)
error->all(FLERR,"Angles assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
void ReadData::dihedrals()
{
int i,m,nchunk;
bigint nread = 0;
bigint ndihedrals = atom->ndihedrals;
while (nread < ndihedrals) {
nchunk = MIN(ndihedrals-nread,CHUNK);
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_dihedrals(nchunk,buffer);
nread += nchunk;
}
// check that dihedrals were assigned correctly
int nlocal = atom->nlocal;
bigint sum;
bigint n = 0;
for (i = 0; i < nlocal; i++) n += atom->num_dihedral[i];
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 4;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " dihedrals\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " dihedrals\n",sum/factor);
}
if (sum != factor*atom->ndihedrals)
error->all(FLERR,"Dihedrals assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
void ReadData::impropers()
{
int i,m,nchunk;
bigint nread = 0;
bigint nimpropers = atom->nimpropers;
while (nread < nimpropers) {
nchunk = MIN(nimpropers-nread,CHUNK);
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_impropers(nchunk,buffer);
nread += nchunk;
}
// check that impropers were assigned correctly
int nlocal = atom->nlocal;
bigint sum;
bigint n = 0;
for (i = 0; i < nlocal; i++) n += atom->num_improper[i];
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 4;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " impropers\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " impropers\n",sum/factor);
}
if (sum != factor*atom->nimpropers)
error->all(FLERR,"Impropers assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
void ReadData::mass()
{
int i,m;
char *buf = new char[atom->ntypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->ntypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->ntypes; i++) {
atom->set_mass(buf);
buf += strlen(buf) + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::paircoeffs()
{
int i,m;
char *buf = new char[atom->ntypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->ntypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->ntypes; i++) {
m = strlen(buf) + 1;
parse_coeffs(buf,NULL,1);
force->pair->coeff(narg,arg);
buf += m;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::bondcoeffs()
{
int i,m;
char *buf = new char[atom->nbondtypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->nbondtypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->nbondtypes; i++) {
m = strlen(buf) + 1;
parse_coeffs(buf,NULL,0);
force->bond->coeff(narg,arg);
buf += m;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::anglecoeffs(int which)
{
int i,m;
char *buf = new char[atom->nangletypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->nangletypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->nangletypes; i++) {
m = strlen(buf) + 1;
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"bb",0);
else if (which == 2) parse_coeffs(buf,"ba",0);
force->angle->coeff(narg,arg);
buf += m;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::dihedralcoeffs(int which)
{
int i,m;
char *buf = new char[atom->ndihedraltypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->ndihedraltypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->ndihedraltypes; i++) {
m = strlen(buf) + 1;
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"mbt",0);
else if (which == 2) parse_coeffs(buf,"ebt",0);
else if (which == 3) parse_coeffs(buf,"at",0);
else if (which == 4) parse_coeffs(buf,"aat",0);
else if (which == 5) parse_coeffs(buf,"bb13",0);
force->dihedral->coeff(narg,arg);
buf += m;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::impropercoeffs(int which)
{
int i,m;
char *buf = new char[atom->nimpropertypes*MAXLINE];
char *original = buf;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < atom->nimpropertypes; i++) {
eof = fgets(&buf[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buf[m]);
- if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
+ if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
buf[m-1] = '\0';
}
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buf,m,MPI_CHAR,0,world);
for (i = 0; i < atom->nimpropertypes; i++) {
m = strlen(buf) + 1;
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"aa",0);
force->improper->coeff(narg,arg);
buf += m;
}
delete [] original;
}
/* ----------------------------------------------------------------------
read fix section, pass lines to fix to process
n = index of fix
------------------------------------------------------------------------- */
void ReadData::fix(int ifix, char *line, bigint nlines)
{
int i,m,nchunk;
bigint nread = 0;
while (nread < nlines) {
if (nlines-nread > CHUNK) nchunk = CHUNK;
else nchunk = nlines-nread;
if (me == 0) {
char *eof;
m = 0;
for (i = 0; i < nchunk; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
modify->fix[ifix]->read_data_section(line,nchunk,buffer);
nread += nchunk;
}
}
/* ----------------------------------------------------------------------
proc 0 scans the data file for topology maximums
------------------------------------------------------------------------- */
void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
int &dihedral_per_atom, int &improper_per_atom)
{
int i,tmp1,tmp2,atom1,atom2,atom3,atom4;
char *eof;
if (atom->natoms > MAXSMALLINT)
error->one(FLERR,"Molecular data file has too many atoms");
// customize for new sections
int natoms = static_cast<int> (atom->natoms);
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
int ellipsoid_flag = 0;
int line_flag = 0;
int tri_flag = 0;
// customize for new sections
// allocate topology counting vector
// initially, array length = 1 to natoms
// will grow via reallocate() if atom IDs > natoms
int cmax = natoms + 1;
int *count;
memory->create(count,cmax,"read_data:count");
while (strlen(keyword)) {
// allow special fixes first chance to match and process the section
// if fix matches, continue to next section
if (nfix) {
for (i = 0; i < nfix; i++) {
printf("LINE SECTION %s %s\n",line,fix_section[i]);
if (strstr(line,fix_section[i])) {
int n = modify->fix[fix_index[i]]->read_data_skip_lines(keyword);
printf("NLINES SKIP %d\n",n);
skip_lines(n);
parse_keyword(0,0);
break;
}
}
if (i < nfix) continue;
}
if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes);
else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms);
else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms);
else if (strcmp(keyword,"Ellipsoids") == 0) {
if (!avec_ellipsoid)
error->one(FLERR,"Invalid data file section: Ellipsoids");
ellipsoid_flag = 1;
skip_lines(nellipsoids);
} else if (strcmp(keyword,"Lines") == 0) {
if (!avec_line) error->one(FLERR,"Invalid data file section: Lines");
line_flag = 1;
skip_lines(nlines);
} else if (strcmp(keyword,"Triangles") == 0) {
if (!avec_tri) error->one(FLERR,"Invalid data file section: Triangles");
tri_flag = 1;
skip_lines(ntris);
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
error->one(FLERR,"Must define pair_style before Pair Coeffs");
skip_lines(atom->ntypes);
} else if (strcmp(keyword,"Bond Coeffs") == 0) {
if (atom->avec->bonds_allow == 0)
error->one(FLERR,"Invalid data file section: Bond Coeffs");
if (force->bond == NULL)
error->one(FLERR,"Must define bond_style before Bond Coeffs");
skip_lines(atom->nbondtypes);
} else if (strcmp(keyword,"Angle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->one(FLERR,"Invalid data file section: Angle Coeffs");
if (force->angle == NULL)
error->one(FLERR,"Must define angle_style before Angle Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"Dihedral Coeffs") == 0) {
skip_lines(atom->ndihedraltypes);
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: Dihedral Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,"Must define dihedral_style before Dihedral Coeffs");
} else if (strcmp(keyword,"Improper Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->one(FLERR,"Invalid data file section: Improper Coeffs");
if (force->improper == NULL)
error->one(FLERR,"Must define improper_style before Improper Coeffs");
skip_lines(atom->nimpropertypes);
} else if (strcmp(keyword,"BondBond Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->one(FLERR,"Invalid data file section: BondBond Coeffs");
if (force->angle == NULL)
error->one(FLERR,"Must define angle_style before BondBond Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"BondAngle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->one(FLERR,"Invalid data file section: BondAngle Coeffs");
if (force->angle == NULL)
error->one(FLERR,"Must define angle_style before BondAngle Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->one(FLERR,"Invalid data file section: BondBond13 Coeffs");
if (force->dihedral == NULL)
error->one(FLERR,"Must define dihedral_style before BondBond13 Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngle Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->one(FLERR,"Invalid data file section: AngleAngle Coeffs");
if (force->improper == NULL)
error->one(FLERR,"Must define improper_style before AngleAngle Coeffs");
skip_lines(atom->nimpropertypes);
} else if (strcmp(keyword,"Bonds") == 0) {
for (i = 1; i < cmax; i++) count[i] = 0;
if (force->newton_bond)
for (i = 0; i < atom->nbonds; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d",&tmp1,&tmp2,&atom1,&atom2);
if (atom1 >= cmax) cmax = reallocate(&count,cmax,atom1);
count[atom1]++;
}
else
for (i = 0; i < atom->nbonds; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d",&tmp1,&tmp2,&atom1,&atom2);
int amax = MAX(atom1,atom2);
if (amax >= cmax) cmax = reallocate(&count,cmax,amax);
count[atom1]++;
count[atom2]++;
}
for (i = 1; i < cmax; i++) bond_per_atom = MAX(bond_per_atom,count[i]);
if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per_atom);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per_atom);
} else if (strcmp(keyword,"Angles") == 0) {
for (i = 1; i < cmax; i++) count[i] = 0;
if (force->newton_bond)
for (i = 0; i < atom->nangles; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d",&tmp1,&tmp2,&atom1,&atom2,&atom3);
if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2);
count[atom2]++;
}
else
for (i = 0; i < atom->nangles; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d",&tmp1,&tmp2,&atom1,&atom2,&atom3);
int amax = MAX(atom1,atom2);
amax = MAX(amax,atom3);
if (amax >= cmax) cmax = reallocate(&count,cmax,amax);
count[atom1]++;
count[atom2]++;
count[atom3]++;
}
for (i = 1; i < cmax; i++) angle_per_atom = MAX(angle_per_atom,count[i]);
if (screen) fprintf(screen," %d = max angles/atom\n",angle_per_atom);
if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per_atom);
} else if (strcmp(keyword,"Dihedrals") == 0) {
for (i = 1; i < cmax; i++) count[i] = 0;
if (force->newton_bond)
for (i = 0; i < atom->ndihedrals; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d %d",
&tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4);
if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2);
count[atom2]++;
}
else
for (i = 0; i < atom->ndihedrals; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d %d",
&tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4);
int amax = MAX(atom1,atom2);
amax = MAX(amax,atom3);
amax = MAX(amax,atom4);
if (amax >= cmax) cmax = reallocate(&count,cmax,amax);
count[atom1]++;
count[atom2]++;
count[atom3]++;
count[atom4]++;
}
for (i = 1; i < cmax; i++)
dihedral_per_atom = MAX(dihedral_per_atom,count[i]);
if (screen)
fprintf(screen," %d = max dihedrals/atom\n",dihedral_per_atom);
if (logfile)
fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per_atom);
} else if (strcmp(keyword,"Impropers") == 0) {
for (i = 1; i < cmax; i++) count[i] = 0;
if (force->newton_bond)
for (i = 0; i < atom->nimpropers; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d %d",
&tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4);
if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2);
count[atom2]++;
}
else
for (i = 0; i < atom->nimpropers; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(line,"%d %d %d %d %d %d",
&tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4);
int amax = MAX(atom1,atom2);
amax = MAX(amax,atom3);
amax = MAX(amax,atom4);
if (amax >= cmax) cmax = reallocate(&count,cmax,amax);
count[atom1]++;
count[atom2]++;
count[atom3]++;
count[atom4]++;
}
for (i = 1; i < cmax; i++)
improper_per_atom = MAX(improper_per_atom,count[i]);
if (screen)
fprintf(screen," %d = max impropers/atom\n",improper_per_atom);
if (logfile)
fprintf(logfile," %d = max impropers/atom\n",improper_per_atom);
} else {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->one(FLERR,str);
}
parse_keyword(0,0);
}
// free topology counting vector
memory->destroy(count);
// error check that topology was specified in file
if ((atom->nbonds && !bond_per_atom) ||
(atom->nangles && !angle_per_atom) ||
(atom->ndihedrals && !dihedral_per_atom) ||
(atom->nimpropers && !improper_per_atom))
error->one(FLERR,"Needed topology not in data file");
// customize for new sections
// error check that Bonus sections were speficied in file
if (nellipsoids && !ellipsoid_flag)
error->one(FLERR,"Needed bonus data not in data file");
if (nlines && !line_flag)
error->one(FLERR,"Needed bonus data not in data file");
if (ntris && !tri_flag)
error->one(FLERR,"Needed bonus data not in data file");
}
/* ----------------------------------------------------------------------
reallocate the count vector from cmax to amax+1 and return new length
zero new locations
------------------------------------------------------------------------- */
int ReadData::reallocate(int **pcount, int cmax, int amax)
{
int *count = *pcount;
memory->grow(count,amax+1,"read_data:count");
for (int i = cmax; i <= amax; i++) count[i] = 0;
*pcount = count;
return amax+1;
}
/* ----------------------------------------------------------------------
proc 0 opens data file
test if gzipped
------------------------------------------------------------------------- */
void ReadData::open(char *file)
{
compressed = 0;
char *suffix = file + strlen(file) - 3;
if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
sprintf(gunzip,"gunzip -c %s",file);
fp = popen(gunzip,"r");
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
}
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
}
/* ----------------------------------------------------------------------
grab next keyword
read lines until one is non-blank
keyword is all text on line w/out leading & trailing white space
read one additional line (assumed blank)
if any read hits EOF, set keyword to empty
if first = 1, line variable holds non-blank line that ended header
if flag = 0, only proc 0 is calling so no bcast
else flag = 1, bcast keyword line to all procs
------------------------------------------------------------------------- */
void ReadData::parse_keyword(int first, int flag)
{
int eof = 0;
// proc 0 reads upto non-blank line plus 1 following line
// eof is set to 1 if any read hits end-of-file
if (me == 0) {
if (!first) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
}
while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
}
if (fgets(buffer,MAXLINE,fp) == NULL) eof = 1;
}
// if eof, set keyword empty and return
if (flag) MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) {
keyword[0] = '\0';
return;
}
// bcast keyword line to all procs
if (flag) {
int n;
if (me == 0) n = strlen(line) + 1;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
}
// copy non-whitespace portion of line into keyword
int start = strspn(line," \t\n\r");
int stop = strlen(line) - 1;
while (line[stop] == ' ' || line[stop] == '\t'
|| line[stop] == '\n' || line[stop] == '\r') stop--;
line[stop+1] = '\0';
strcpy(keyword,&line[start]);
}
/* ----------------------------------------------------------------------
proc 0 reads N lines from file
------------------------------------------------------------------------- */
void ReadData::skip_lines(int n)
{
char *eof;
for (int i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
}
/* ----------------------------------------------------------------------
parse a line of coeffs into words, storing them in narg,arg
trim anything from '#' onward
word strings remain in line, are not copied
if addstr != NULL, add addstr as extra arg for class2 angle/dihedral/improper
if 2nd word starts with letter, then is hybrid style, add addstr after it
else add addstr before 2nd word
if dupflag, duplicate 1st word, so pair_coeff "2" becomes "2 2"
------------------------------------------------------------------------- */
void ReadData::parse_coeffs(char *line, const char *addstr, int dupflag)
{
char *ptr;
if (ptr = strchr(line,'#')) *ptr = '\0';
narg = 0;
char *word = strtok(line," \t\n\r\f");
while (word) {
if (narg == maxarg) {
maxarg += DELTA;
arg = (char **)
memory->srealloc(arg,maxarg*sizeof(char *),"read_data:arg");
}
if (addstr && narg == 1 && !islower(word[0])) arg[narg++] = (char *) addstr;
arg[narg++] = word;
if (addstr && narg == 2 && islower(word[0])) arg[narg++] = (char *) addstr;
if (dupflag && narg == 1) arg[narg++] = word;
word = strtok(NULL," \t\n\r\f");
}
}