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

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

This commit is contained in:
sjplimp 2013-06-27 23:50:33 +00:00
parent 1e4b3bcf87
commit 440d389779
4 changed files with 4378 additions and 4367 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6784
src/KSPACE/pppm.cpp
View File

File diff suppressed because it is too large Load Diff

670
src/KSPACE/pppm.h
View File

@ -1,335 +1,335 @@
/* -*- 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 h_x,h_y,h_z;
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;
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();
virtual void compute_gf_ik();
virtual void compute_gf_ad();
void compute_sf_precoeff();
virtual void particle_map();
virtual void make_rho();
virtual void brick2fft();
virtual void poisson();
virtual void poisson_ik();
virtual void poisson_ad();
virtual void fieldforce();
virtual void fieldforce_ik();
virtual 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
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 *);
// triclinic
int triclinic; // domain settings, orthog or triclinic
void setup_triclinic();
void compute_gf_ik_triclinic();
void poisson_ik_triclinic();
void poisson_groups_triclinic();
// 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 and 'kspace_modify diff ad'
This feature is not yet supported.
E: Cannot (yet) use PPPM with triclinic box and slab correction
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 or dispersion component be used.
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 (yet) use PPPM with triclinic box and TIP4P
This feature is not yet supported.
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 the long-range Coulombic solvers.
W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor
This may lead to a larger grid than desired. See the kspace_modify overlap
command to prevent changing of the PPPM order.
E: PPPM order < minimum allowed order
The default minimum order is 2. This can be reset by the
kspace_modify minorder command.
E: PPPM grid stencil extends beyond nearest neighbor processor
This is not allowed if the kspace_modify overlap setting is no.
E: KSpace accuracy must be > 0
The kspace accuracy designated in the input must be greater than zero.
E: Could not compute grid size
The code is unable to compute a grid size consistent with the desired
accuracy. This error should not occur for typical problems. Please
send an email to the developers.
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: Could not compute g_ewald
The Newton-Raphson solver failed to converge to a good value for
g_ewald. This error should not occur for typical problems. Please
send an email to the developers.
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.
E: Cannot (yet) use 'kspace_modify diff ad' with compute group/group
This option is not yet supported.
*/
/* -*- 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 h_x,h_y,h_z;
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;
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);
virtual double compute_qopt();
virtual void compute_gf_denom();
virtual void compute_gf_ik();
virtual void compute_gf_ad();
void compute_sf_precoeff();
virtual void particle_map();
virtual void make_rho();
virtual void brick2fft();
virtual void poisson();
virtual void poisson_ik();
virtual void poisson_ad();
virtual void fieldforce();
virtual void fieldforce_ik();
virtual 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
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 *);
// triclinic
int triclinic; // domain settings, orthog or triclinic
void setup_triclinic();
void compute_gf_ik_triclinic();
void poisson_ik_triclinic();
void poisson_groups_triclinic();
// 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 and 'kspace_modify diff ad'
This feature is not yet supported.
E: Cannot (yet) use PPPM with triclinic box and slab correction
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 or dispersion component be used.
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 (yet) use PPPM with triclinic box and TIP4P
This feature is not yet supported.
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 the long-range Coulombic solvers.
W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor
This may lead to a larger grid than desired. See the kspace_modify overlap
command to prevent changing of the PPPM order.
E: PPPM order < minimum allowed order
The default minimum order is 2. This can be reset by the
kspace_modify minorder command.
E: PPPM grid stencil extends beyond nearest neighbor processor
This is not allowed if the kspace_modify overlap setting is no.
E: KSpace accuracy must be > 0
The kspace accuracy designated in the input must be greater than zero.
E: Could not compute grid size
The code is unable to compute a grid size consistent with the desired
accuracy. This error should not occur for typical problems. Please
send an email to the developers.
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: Could not compute g_ewald
The Newton-Raphson solver failed to converge to a good value for
g_ewald. This error should not occur for typical problems. Please
send an email to the developers.
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.
E: Cannot (yet) use 'kspace_modify diff ad' with compute group/group
This option is not yet supported.
*/

893
src/kspace.cpp
View File

@ -1,446 +1,447 @@
/* ----------------------------------------------------------------------
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 "stdlib.h"
#include "string.h"
#include "kspace.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "memory.h"
#include "atom_masks.h"
#include "error.h"
#include "suffix.h"
#include "domain.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
{
energy = 0.0;
virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
triclinic_support = 1;
ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = 0;
compute_flag = 1;
group_group_enable = 0;
order = 5;
gridflag = 0;
gewaldflag = 0;
minorder = 2;
overlap_allowed = 1;
fftbench = 1;
order_6 = 5;
gridflag_6 = 0;
gewaldflag_6 = 0;
slabflag = 0;
differentiation_flag = 0;
slab_volfactor = 1;
suffix_flag = Suffix::NONE;
adjust_cutoff_flag = 1;
accuracy_absolute = -1.0;
two_charge_force = force->qqr2e *
(force->qelectron * force->qelectron) /
(force->angstrom * force->angstrom);
maxeatom = maxvatom = 0;
eatom = NULL;
vatom = NULL;
datamask = ALL_MASK;
datamask_ext = ALL_MASK;
memory->create(gcons,7,7,"kspace:gcons");
gcons[2][0] = 15.0 / 8.0;
gcons[2][1] = -5.0 / 4.0;
gcons[2][2] = 3.0 / 8.0;
gcons[3][0] = 35.0 / 16.0;
gcons[3][1] = -35.0 / 16.0;
gcons[3][2] = 21.0 / 16.0;
gcons[3][3] = -5.0 / 16.0;
gcons[4][0] = 315.0 / 128.0;
gcons[4][1] = -105.0 / 32.0;
gcons[4][2] = 189.0 / 64.0;
gcons[4][3] = -45.0 / 32.0;
gcons[4][4] = 35.0 / 128.0;
gcons[5][0] = 693.0 / 256.0;
gcons[5][1] = -1155.0 / 256.0;
gcons[5][2] = 693.0 / 128.0;
gcons[5][3] = -495.0 / 128.0;
gcons[5][4] = 385.0 / 256.0;
gcons[5][5] = -63.0 / 256.0;
gcons[6][0] = 3003.0 / 1024.0;
gcons[6][1] = -3003.0 / 512.0;
gcons[6][2] = 9009.0 / 1024.0;
gcons[6][3] = -2145.0 / 256.0;
gcons[6][4] = 5005.0 / 1024.0;
gcons[6][5] = -819.0 / 512.0;
gcons[6][6] = 231.0 / 1024.0;
memory->create(dgcons,7,6,"kspace:dgcons");
dgcons[2][0] = -5.0 / 2.0;
dgcons[2][1] = 3.0 / 2.0;
dgcons[3][0] = -35.0 / 8.0;
dgcons[3][1] = 21.0 / 4.0;
dgcons[3][2] = -15.0 / 8.0;
dgcons[4][0] = -105.0 / 16.0;
dgcons[4][1] = 189.0 / 16.0;
dgcons[4][2] = -135.0 / 16.0;
dgcons[4][3] = 35.0 / 16.0;
dgcons[5][0] = -1155.0 / 128.0;
dgcons[5][1] = 693.0 / 32.0;
dgcons[5][2] = -1485.0 / 64.0;
dgcons[5][3] = 385.0 / 32.0;
dgcons[5][4] = -315.0 / 128.0;
dgcons[6][0] = -3003.0 / 256.0;
dgcons[6][1] = 9009.0 / 256.0;
dgcons[6][2] = -6435.0 / 128.0;
dgcons[6][3] = 5005.0 / 128.0;
dgcons[6][4] = -4095.0 / 256.0;
dgcons[6][5] = 693.0 / 256.0;
}
/* ---------------------------------------------------------------------- */
KSpace::~KSpace()
{
memory->destroy(eatom);
memory->destroy(vatom);
memory->destroy(gcons);
memory->destroy(dgcons);
}
/* ---------------------------------------------------------------------- */
void KSpace::triclinic_check()
{
if (domain->triclinic && triclinic_support != 1)
error->all(FLERR,"KSpace style does not yet support triclinic geometries");
}
/* ---------------------------------------------------------------------- */
void KSpace::compute_dummy(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
}
/* ----------------------------------------------------------------------
check that pair style is compatible with long-range solver
------------------------------------------------------------------------- */
void KSpace::pair_check()
{
if (force->pair == NULL)
error->all(FLERR,"KSpace solver requires a pair style");
if (ewaldflag && force->pair->ewaldflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (pppmflag && force->pair->pppmflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (msmflag && force->pair->msmflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (dispersionflag && force->pair->dispersionflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (tip4pflag && force->pair->tip4pflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (dipoleflag && force->pair->dipoleflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
}
/* ----------------------------------------------------------------------
setup for energy, virial computation
see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
------------------------------------------------------------------------- */
void KSpace::ev_setup(int eflag, int vflag)
{
int i,n;
evflag = 1;
eflag_either = eflag;
eflag_global = eflag % 2;
eflag_atom = eflag / 2;
vflag_either = vflag;
vflag_global = vflag % 4;
vflag_atom = vflag / 4;
if (eflag_atom || vflag_atom) evflag_atom = 1;
else evflag_atom = 0;
// reallocate per-atom arrays if necessary
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
memory->create(eatom,maxeatom,"kspace:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
memory->create(vatom,maxvatom,6,"kspace:vatom");
}
// zero accumulators
if (eflag_global) energy = 0.0;
if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
if (eflag_atom) {
n = atom->nlocal;
if (tip4pflag) n += atom->nghost;
for (i = 0; i < n; i++) eatom[i] = 0.0;
}
if (vflag_atom) {
n = atom->nlocal;
if (tip4pflag) n += atom->nghost;
for (i = 0; i < n; i++) {
vatom[i][0] = 0.0;
vatom[i][1] = 0.0;
vatom[i][2] = 0.0;
vatom[i][3] = 0.0;
vatom[i][4] = 0.0;
vatom[i][5] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
estimate the accuracy of the short-range coulomb tables
------------------------------------------------------------------------- */
double KSpace::estimate_table_accuracy(double q2_over_sqrt, double spr)
{
double table_accuracy = 0.0;
int nctb = force->pair->ncoultablebits;
if (nctb) {
double empirical_precision[17];
empirical_precision[6] = 6.99E-03;
empirical_precision[7] = 1.78E-03;
empirical_precision[8] = 4.72E-04;
empirical_precision[9] = 1.17E-04;
empirical_precision[10] = 2.95E-05;
empirical_precision[11] = 7.41E-06;
empirical_precision[12] = 1.76E-06;
empirical_precision[13] = 9.28E-07;
empirical_precision[14] = 7.46E-07;
empirical_precision[15] = 7.32E-07;
empirical_precision[16] = 7.30E-07;
if (nctb <= 6) table_accuracy = empirical_precision[6];
else if (nctb <= 16) table_accuracy = empirical_precision[nctb];
else table_accuracy = empirical_precision[16];
table_accuracy *= q2_over_sqrt;
if ((table_accuracy > spr) && (comm->me == 0))
error->warning(FLERR,"For better accuracy use 'pair_modify table 0'");
}
return table_accuracy;
}
/* ----------------------------------------------------------------------
convert box coords vector to transposed triclinic lamda (0-1) coords
vector, lamda = [(H^-1)^T] * v, does not preserve vector magnitude
v and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::x2lamdaT(double *v, double *lamda)
{
double *h_inv = domain->h_inv;
double lamda_tmp[3];
lamda_tmp[0] = h_inv[0]*v[0];
lamda_tmp[1] = h_inv[5]*v[0] + h_inv[1]*v[1];
lamda_tmp[2] = h_inv[4]*v[0] + h_inv[3]*v[1] + h_inv[2]*v[2];
lamda[0] = lamda_tmp[0];
lamda[1] = lamda_tmp[1];
lamda[2] = lamda_tmp[2];
}
/* ----------------------------------------------------------------------
convert lamda (0-1) coords vector to transposed box coords vector
lamda = (H^T) * v, does not preserve vector magnitude
v and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::lamda2xT(double *lamda, double *v)
{
double *h = domain->h;
double v_tmp[3];
v_tmp[0] = h[0]*lamda[0];
v_tmp[1] = h[5]*lamda[0] + h[1]*lamda[1];
v_tmp[2] = h[4]*lamda[0] + h[3]*lamda[1] + h[2]*lamda[2];
v[0] = v_tmp[0];
v[1] = v_tmp[1];
v[2] = v_tmp[2];
}
/* ----------------------------------------------------------------------
convert triclinic lamda (0-1) coords vector to box coords vector
v = H * lamda, does not preserve vector magnitude
lamda and v can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::lamda2xvector(double *lamda, double *v)
{
double *h = domain->h;
v[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2];
v[1] = h[1]*lamda[1] + h[3]*lamda[2];
v[2] = h[2]*lamda[2];
}
/* ----------------------------------------------------------------------
convert a sphere in box coords to an ellipsoid in lamda (0-1)
coords and return the tight (axis-aligned) bounding box, does not
preserve vector magnitude
see http://www.loria.fr/~shornus/ellipsoid-bbox.html and
http://yiningkarlli.blogspot.com/2013/02/bounding-boxes-for-ellipsoidsfigure.html
------------------------------------------------------------------------- */
void KSpace::kspacebbox(double r, double *b)
{
double *h = domain->h;
double lx,ly,lz,xy,xz,yz;
lx = h[0];
ly = h[1];
lz = h[2];
yz = h[3];
xz = h[4];
xy = h[5];
b[0] = r*sqrt(ly*ly*lz*lz + ly*ly*xz*xz - 2.0*ly*xy*xz*yz + lz*lz*xy*xy +
xy*xy*yz*yz)/(lx*ly*lz);
b[1] = r*sqrt(lz*lz + yz*yz)/(ly*lz);
b[2] = r/lz;
}
/* ----------------------------------------------------------------------
modify parameters of the KSpace style
------------------------------------------------------------------------- */
void KSpace::modify_params(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"mesh") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
nx_pppm = nx_msm_max = atoi(arg[iarg+1]);
ny_pppm = ny_msm_max = atoi(arg[iarg+2]);
nz_pppm = nz_msm_max = atoi(arg[iarg+3]);
if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0) gridflag = 0;
else gridflag = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"mesh/disp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
nx_pppm_6 = atoi(arg[iarg+1]);
ny_pppm_6 = atoi(arg[iarg+2]);
nz_pppm_6 = atoi(arg[iarg+3]);
if (nx_pppm_6 == 0 || ny_pppm_6 == 0 || nz_pppm_6 == 0) gridflag_6 = 0;
else gridflag_6 = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"order") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"order/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order_6 = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"minorder") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
minorder = atoi(arg[iarg+1]);
if (minorder < 2) error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"overlap") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) overlap_allowed = 1;
else if (strcmp(arg[iarg+1],"no") == 0) overlap_allowed = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"force") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
accuracy_absolute = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"gewald") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
g_ewald = atof(arg[iarg+1]);
if (g_ewald == 0.0) gewaldflag = 0;
else gewaldflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"gewald/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
g_ewald_6 = atof(arg[iarg+1]);
if (g_ewald_6 == 0.0) gewaldflag_6 = 0;
else gewaldflag_6 = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"slab") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"nozforce") == 0) {
slabflag = 2;
} else {
slabflag = 1;
slab_volfactor = atof(arg[iarg+1]);
if (slab_volfactor <= 1.0)
error->all(FLERR,"Bad kspace_modify slab parameter");
if (slab_volfactor < 2.0 && comm->me == 0)
error->warning(FLERR,"Kspace_modify slab param < 2.0 may "
"cause unphysical behavior");
}
iarg += 2;
} else if (strcmp(arg[iarg],"compute") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"fftbench") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) fftbench = 1;
else if (strcmp(arg[iarg+1],"no") == 0) fftbench = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"diff") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"ad") == 0) differentiation_flag = 1;
else if (strcmp(arg[iarg+1],"ik") == 0) differentiation_flag = 0;
else error->all(FLERR, "Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"cutoff/adjust") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) adjust_cutoff_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) adjust_cutoff_flag = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else error->all(FLERR,"Illegal kspace_modify command");
}
}
/* ---------------------------------------------------------------------- */
void *KSpace::extract(const char *str)
{
if (strcmp(str,"scale") == 0) return (void *) &scale;
return NULL;
}
/* ----------------------------------------------------------------------
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 "stdlib.h"
#include "string.h"
#include "kspace.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "memory.h"
#include "atom_masks.h"
#include "error.h"
#include "suffix.h"
#include "domain.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
{
energy = 0.0;
virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
triclinic_support = 1;
ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = 0;
compute_flag = 1;
group_group_enable = 0;
stagger_flag = 0;
order = 5;
gridflag = 0;
gewaldflag = 0;
minorder = 2;
overlap_allowed = 1;
fftbench = 1;
order_6 = 5;
gridflag_6 = 0;
gewaldflag_6 = 0;
slabflag = 0;
differentiation_flag = 0;
slab_volfactor = 1;
suffix_flag = Suffix::NONE;
adjust_cutoff_flag = 1;
accuracy_absolute = -1.0;
two_charge_force = force->qqr2e *
(force->qelectron * force->qelectron) /
(force->angstrom * force->angstrom);
maxeatom = maxvatom = 0;
eatom = NULL;
vatom = NULL;
datamask = ALL_MASK;
datamask_ext = ALL_MASK;
memory->create(gcons,7,7,"kspace:gcons");
gcons[2][0] = 15.0 / 8.0;
gcons[2][1] = -5.0 / 4.0;
gcons[2][2] = 3.0 / 8.0;
gcons[3][0] = 35.0 / 16.0;
gcons[3][1] = -35.0 / 16.0;
gcons[3][2] = 21.0 / 16.0;
gcons[3][3] = -5.0 / 16.0;
gcons[4][0] = 315.0 / 128.0;
gcons[4][1] = -105.0 / 32.0;
gcons[4][2] = 189.0 / 64.0;
gcons[4][3] = -45.0 / 32.0;
gcons[4][4] = 35.0 / 128.0;
gcons[5][0] = 693.0 / 256.0;
gcons[5][1] = -1155.0 / 256.0;
gcons[5][2] = 693.0 / 128.0;
gcons[5][3] = -495.0 / 128.0;
gcons[5][4] = 385.0 / 256.0;
gcons[5][5] = -63.0 / 256.0;
gcons[6][0] = 3003.0 / 1024.0;
gcons[6][1] = -3003.0 / 512.0;
gcons[6][2] = 9009.0 / 1024.0;
gcons[6][3] = -2145.0 / 256.0;
gcons[6][4] = 5005.0 / 1024.0;
gcons[6][5] = -819.0 / 512.0;
gcons[6][6] = 231.0 / 1024.0;
memory->create(dgcons,7,6,"kspace:dgcons");
dgcons[2][0] = -5.0 / 2.0;
dgcons[2][1] = 3.0 / 2.0;
dgcons[3][0] = -35.0 / 8.0;
dgcons[3][1] = 21.0 / 4.0;
dgcons[3][2] = -15.0 / 8.0;
dgcons[4][0] = -105.0 / 16.0;
dgcons[4][1] = 189.0 / 16.0;
dgcons[4][2] = -135.0 / 16.0;
dgcons[4][3] = 35.0 / 16.0;
dgcons[5][0] = -1155.0 / 128.0;
dgcons[5][1] = 693.0 / 32.0;
dgcons[5][2] = -1485.0 / 64.0;
dgcons[5][3] = 385.0 / 32.0;
dgcons[5][4] = -315.0 / 128.0;
dgcons[6][0] = -3003.0 / 256.0;
dgcons[6][1] = 9009.0 / 256.0;
dgcons[6][2] = -6435.0 / 128.0;
dgcons[6][3] = 5005.0 / 128.0;
dgcons[6][4] = -4095.0 / 256.0;
dgcons[6][5] = 693.0 / 256.0;
}
/* ---------------------------------------------------------------------- */
KSpace::~KSpace()
{
memory->destroy(eatom);
memory->destroy(vatom);
memory->destroy(gcons);
memory->destroy(dgcons);
}
/* ---------------------------------------------------------------------- */
void KSpace::triclinic_check()
{
if (domain->triclinic && triclinic_support != 1)
error->all(FLERR,"KSpace style does not yet support triclinic geometries");
}
/* ---------------------------------------------------------------------- */
void KSpace::compute_dummy(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
}
/* ----------------------------------------------------------------------
check that pair style is compatible with long-range solver
------------------------------------------------------------------------- */
void KSpace::pair_check()
{
if (force->pair == NULL)
error->all(FLERR,"KSpace solver requires a pair style");
if (ewaldflag && force->pair->ewaldflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (pppmflag && force->pair->pppmflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (msmflag && force->pair->msmflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (dispersionflag && force->pair->dispersionflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (tip4pflag && force->pair->tip4pflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
if (dipoleflag && force->pair->dipoleflag == 0)
error->all(FLERR,"KSpace style is incompatible with Pair style");
}
/* ----------------------------------------------------------------------
setup for energy, virial computation
see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
------------------------------------------------------------------------- */
void KSpace::ev_setup(int eflag, int vflag)
{
int i,n;
evflag = 1;
eflag_either = eflag;
eflag_global = eflag % 2;
eflag_atom = eflag / 2;
vflag_either = vflag;
vflag_global = vflag % 4;
vflag_atom = vflag / 4;
if (eflag_atom || vflag_atom) evflag_atom = 1;
else evflag_atom = 0;
// reallocate per-atom arrays if necessary
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
memory->create(eatom,maxeatom,"kspace:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
memory->create(vatom,maxvatom,6,"kspace:vatom");
}
// zero accumulators
if (eflag_global) energy = 0.0;
if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
if (eflag_atom) {
n = atom->nlocal;
if (tip4pflag) n += atom->nghost;
for (i = 0; i < n; i++) eatom[i] = 0.0;
}
if (vflag_atom) {
n = atom->nlocal;
if (tip4pflag) n += atom->nghost;
for (i = 0; i < n; i++) {
vatom[i][0] = 0.0;
vatom[i][1] = 0.0;
vatom[i][2] = 0.0;
vatom[i][3] = 0.0;
vatom[i][4] = 0.0;
vatom[i][5] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
estimate the accuracy of the short-range coulomb tables
------------------------------------------------------------------------- */
double KSpace::estimate_table_accuracy(double q2_over_sqrt, double spr)
{
double table_accuracy = 0.0;
int nctb = force->pair->ncoultablebits;
if (nctb) {
double empirical_precision[17];
empirical_precision[6] = 6.99E-03;
empirical_precision[7] = 1.78E-03;
empirical_precision[8] = 4.72E-04;
empirical_precision[9] = 1.17E-04;
empirical_precision[10] = 2.95E-05;
empirical_precision[11] = 7.41E-06;
empirical_precision[12] = 1.76E-06;
empirical_precision[13] = 9.28E-07;
empirical_precision[14] = 7.46E-07;
empirical_precision[15] = 7.32E-07;
empirical_precision[16] = 7.30E-07;
if (nctb <= 6) table_accuracy = empirical_precision[6];
else if (nctb <= 16) table_accuracy = empirical_precision[nctb];
else table_accuracy = empirical_precision[16];
table_accuracy *= q2_over_sqrt;
if ((table_accuracy > spr) && (comm->me == 0))
error->warning(FLERR,"For better accuracy use 'pair_modify table 0'");
}
return table_accuracy;
}
/* ----------------------------------------------------------------------
convert box coords vector to transposed triclinic lamda (0-1) coords
vector, lamda = [(H^-1)^T] * v, does not preserve vector magnitude
v and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::x2lamdaT(double *v, double *lamda)
{
double *h_inv = domain->h_inv;
double lamda_tmp[3];
lamda_tmp[0] = h_inv[0]*v[0];
lamda_tmp[1] = h_inv[5]*v[0] + h_inv[1]*v[1];
lamda_tmp[2] = h_inv[4]*v[0] + h_inv[3]*v[1] + h_inv[2]*v[2];
lamda[0] = lamda_tmp[0];
lamda[1] = lamda_tmp[1];
lamda[2] = lamda_tmp[2];
}
/* ----------------------------------------------------------------------
convert lamda (0-1) coords vector to transposed box coords vector
lamda = (H^T) * v, does not preserve vector magnitude
v and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::lamda2xT(double *lamda, double *v)
{
double *h = domain->h;
double v_tmp[3];
v_tmp[0] = h[0]*lamda[0];
v_tmp[1] = h[5]*lamda[0] + h[1]*lamda[1];
v_tmp[2] = h[4]*lamda[0] + h[3]*lamda[1] + h[2]*lamda[2];
v[0] = v_tmp[0];
v[1] = v_tmp[1];
v[2] = v_tmp[2];
}
/* ----------------------------------------------------------------------
convert triclinic lamda (0-1) coords vector to box coords vector
v = H * lamda, does not preserve vector magnitude
lamda and v can point to same 3-vector
------------------------------------------------------------------------- */
void KSpace::lamda2xvector(double *lamda, double *v)
{
double *h = domain->h;
v[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2];
v[1] = h[1]*lamda[1] + h[3]*lamda[2];
v[2] = h[2]*lamda[2];
}
/* ----------------------------------------------------------------------
convert a sphere in box coords to an ellipsoid in lamda (0-1)
coords and return the tight (axis-aligned) bounding box, does not
preserve vector magnitude
see http://www.loria.fr/~shornus/ellipsoid-bbox.html and
http://yiningkarlli.blogspot.com/2013/02/bounding-boxes-for-ellipsoidsfigure.html
------------------------------------------------------------------------- */
void KSpace::kspacebbox(double r, double *b)
{
double *h = domain->h;
double lx,ly,lz,xy,xz,yz;
lx = h[0];
ly = h[1];
lz = h[2];
yz = h[3];
xz = h[4];
xy = h[5];
b[0] = r*sqrt(ly*ly*lz*lz + ly*ly*xz*xz - 2.0*ly*xy*xz*yz + lz*lz*xy*xy +
xy*xy*yz*yz)/(lx*ly*lz);
b[1] = r*sqrt(lz*lz + yz*yz)/(ly*lz);
b[2] = r/lz;
}
/* ----------------------------------------------------------------------
modify parameters of the KSpace style
------------------------------------------------------------------------- */
void KSpace::modify_params(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"mesh") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
nx_pppm = nx_msm_max = atoi(arg[iarg+1]);
ny_pppm = ny_msm_max = atoi(arg[iarg+2]);
nz_pppm = nz_msm_max = atoi(arg[iarg+3]);
if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0) gridflag = 0;
else gridflag = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"mesh/disp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
nx_pppm_6 = atoi(arg[iarg+1]);
ny_pppm_6 = atoi(arg[iarg+2]);
nz_pppm_6 = atoi(arg[iarg+3]);
if (nx_pppm_6 == 0 || ny_pppm_6 == 0 || nz_pppm_6 == 0) gridflag_6 = 0;
else gridflag_6 = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"order") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"order/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order_6 = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"minorder") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
minorder = atoi(arg[iarg+1]);
if (minorder < 2) error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"overlap") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) overlap_allowed = 1;
else if (strcmp(arg[iarg+1],"no") == 0) overlap_allowed = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"force") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
accuracy_absolute = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"gewald") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
g_ewald = atof(arg[iarg+1]);
if (g_ewald == 0.0) gewaldflag = 0;
else gewaldflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"gewald/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
g_ewald_6 = atof(arg[iarg+1]);
if (g_ewald_6 == 0.0) gewaldflag_6 = 0;
else gewaldflag_6 = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"slab") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"nozforce") == 0) {
slabflag = 2;
} else {
slabflag = 1;
slab_volfactor = atof(arg[iarg+1]);
if (slab_volfactor <= 1.0)
error->all(FLERR,"Bad kspace_modify slab parameter");
if (slab_volfactor < 2.0 && comm->me == 0)
error->warning(FLERR,"Kspace_modify slab param < 2.0 may "
"cause unphysical behavior");
}
iarg += 2;
} else if (strcmp(arg[iarg],"compute") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"fftbench") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) fftbench = 1;
else if (strcmp(arg[iarg+1],"no") == 0) fftbench = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"diff") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"ad") == 0) differentiation_flag = 1;
else if (strcmp(arg[iarg+1],"ik") == 0) differentiation_flag = 0;
else error->all(FLERR, "Illegal kspace_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"cutoff/adjust") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) adjust_cutoff_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) adjust_cutoff_flag = 0;
else error->all(FLERR,"Illegal kspace_modify command");
iarg += 2;
} else error->all(FLERR,"Illegal kspace_modify command");
}
}
/* ---------------------------------------------------------------------- */
void *KSpace::extract(const char *str)
{
if (strcmp(str,"scale") == 0) return (void *) &scale;
return NULL;
}

398
src/kspace.h
View File

@ -1,198 +1,200 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifndef LMP_KSPACE_H
#define LMP_KSPACE_H
#include "pointers.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
namespace LAMMPS_NS {
class KSpace : protected Pointers {
friend class ThrOMP;
friend class FixOMP;
public:
double energy; // accumulated energies
double energy_1,energy_6;
double virial[6]; // accumlated virial
double *eatom,**vatom; // accumulated per-atom energy/virial
double e2group; // accumulated group-group energy
double f2group[3]; // accumulated group-group force
int triclinic_support; // 1 if supports triclinic geometries
int ewaldflag; // 1 if a Ewald solver
int pppmflag; // 1 if a PPPM solver
int msmflag; // 1 if a MSM solver
int dispersionflag; // 1 if a LJ/dispersion solver
int tip4pflag; // 1 if a TIP4P solver
int dipoleflag; // 1 if a dipole solver
int differentiation_flag;
int slabflag;
double slab_volfactor;
int order,order_6;
double accuracy; // accuracy of KSpace solver (force units)
double accuracy_absolute; // user-specifed accuracy in force units
double accuracy_relative; // user-specified dimensionless accuracy
// accurary = acc_rel * two_charge_force
double two_charge_force; // force in user units of two point
// charges separated by 1 Angstrom
double g_ewald,g_ewald_6;
int nx_pppm,ny_pppm,nz_pppm; // global FFT grid for Coulombics
int nx_pppm_6,ny_pppm_6,nz_pppm_6; // global FFT grid for dispersion
int nx_msm_max,ny_msm_max,nz_msm_max;
int group_group_enable; // 1 if style supports group/group calculation
unsigned int datamask;
unsigned int datamask_ext;
int compute_flag; // 0 if skip compute()
int fftbench; // 0 if skip FFT timing
KSpace(class LAMMPS *, int, char **);
virtual ~KSpace();
void triclinic_check();
void modify_params(int, char **);
void *extract(const char *);
void compute_dummy(int, int);
// triclinic
void x2lamdaT(double *, double *);
void lamda2xT(double *, double *);
void lamda2xvector(double *, double *);
void kspacebbox(double, double *);
// general child-class methods
virtual void init() = 0;
virtual void setup() = 0;
virtual void setup_grid() {};
virtual void compute(int, int) = 0;
virtual void compute_group_group(int, 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 *) {};
virtual int timing(int, double &, double &) {return 0;}
virtual int timing_1d(int, double &) {return 0;}
virtual int timing_3d(int, double &) {return 0;}
virtual double memory_usage() {return 0.0;}
/* ----------------------------------------------------------------------
compute gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164–177
------------------------------------------------------------------------- */
double gamma(const double &rho) const {
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double g = gcons[split_order][0];
double rho_n = rho2;
for (int n=1; n<=split_order; n++) {
g += gcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return g;
} else
return (1.0/rho);
}
/* ----------------------------------------------------------------------
compute the derivative of gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164-177
------------------------------------------------------------------------- */
double dgamma(const double &rho) const {
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double dg = dgcons[split_order][0]*rho;
double rho_n = rho*rho2;
for (int n=1; n<split_order; n++) {
dg += dgcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return dg;
} else
return (-1.0/rho/rho);
}
protected:
int gridflag,gridflag_6;
int gewaldflag,gewaldflag_6;
int minorder,overlap_allowed;
int adjust_cutoff_flag;
int suffix_flag; // suffix compatibility flag
double scale;
double **gcons,**dgcons; // accumulated per-atom energy/virial
int evflag,evflag_atom;
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int maxeatom,maxvatom;
void pair_check();
void ev_setup(int, int);
double estimate_table_accuracy(double, double);
};
}
#endif
/* ERROR/WARNING messages:
E: KSpace solver requires a pair style
No pair style is defined.
E: KSpace style is incompatible with Pair style
Setting a kspace style requires that a pair style with a long-range
Coulombic or dispersion component be used.
W: For better accuracy use 'pair_modify table 0'
The user-specified force accuracy cannot be achieved unless the table
feature is disabled by using 'pair_modify table 0'.
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: Bad kspace_modify slab parameter
Kspace_modify value for the slab/volume keyword must be >= 2.0.
W: Kspace_modify slab param < 2.0 may cause unphysical behavior
The kspace_modify slab parameter should be larger to insure periodic
grids padded with empty space do not overlap.
*/
/* -*- 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.
------------------------------------------------------------------------- */
#ifndef LMP_KSPACE_H
#define LMP_KSPACE_H
#include "pointers.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
namespace LAMMPS_NS {
class KSpace : protected Pointers {
friend class ThrOMP;
friend class FixOMP;
public:
double energy; // accumulated energies
double energy_1,energy_6;
double virial[6]; // accumlated virial
double *eatom,**vatom; // accumulated per-atom energy/virial
double e2group; // accumulated group-group energy
double f2group[3]; // accumulated group-group force
int triclinic_support; // 1 if supports triclinic geometries
int ewaldflag; // 1 if a Ewald solver
int pppmflag; // 1 if a PPPM solver
int msmflag; // 1 if a MSM solver
int dispersionflag; // 1 if a LJ/dispersion solver
int tip4pflag; // 1 if a TIP4P solver
int dipoleflag; // 1 if a dipole solver
int differentiation_flag;
int slabflag;
double slab_volfactor;
int order,order_6;
double accuracy; // accuracy of KSpace solver (force units)
double accuracy_absolute; // user-specifed accuracy in force units
double accuracy_relative; // user-specified dimensionless accuracy
// accurary = acc_rel * two_charge_force
double two_charge_force; // force in user units of two point
// charges separated by 1 Angstrom
double g_ewald,g_ewald_6;
int nx_pppm,ny_pppm,nz_pppm; // global FFT grid for Coulombics
int nx_pppm_6,ny_pppm_6,nz_pppm_6; // global FFT grid for dispersion
int nx_msm_max,ny_msm_max,nz_msm_max;
int group_group_enable; // 1 if style supports group/group calculation
unsigned int datamask;
unsigned int datamask_ext;
int compute_flag; // 0 if skip compute()
int fftbench; // 0 if skip FFT timing
int stagger_flag; // 1 if using staggered PPPM grids
KSpace(class LAMMPS *, int, char **);
virtual ~KSpace();
void triclinic_check();
void modify_params(int, char **);
void *extract(const char *);
void compute_dummy(int, int);
// triclinic
void x2lamdaT(double *, double *);
void lamda2xT(double *, double *);
void lamda2xvector(double *, double *);
void kspacebbox(double, double *);
// general child-class methods
virtual void init() = 0;
virtual void setup() = 0;
virtual void setup_grid() {};
virtual void compute(int, int) = 0;
virtual void compute_group_group(int, 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 *) {};
virtual int timing(int, double &, double &) {return 0;}
virtual int timing_1d(int, double &) {return 0;}
virtual int timing_3d(int, double &) {return 0;}
virtual double memory_usage() {return 0.0;}
/* ----------------------------------------------------------------------
compute gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164–177
------------------------------------------------------------------------- */
double gamma(const double &rho) const {
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double g = gcons[split_order][0];
double rho_n = rho2;
for (int n=1; n<=split_order; n++) {
g += gcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return g;
} else
return (1.0/rho);
}
/* ----------------------------------------------------------------------
compute the derivative of gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164-177
------------------------------------------------------------------------- */
double dgamma(const double &rho) const {
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double dg = dgcons[split_order][0]*rho;
double rho_n = rho*rho2;
for (int n=1; n<split_order; n++) {
dg += dgcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return dg;
} else
return (-1.0/rho/rho);
}
protected:
int gridflag,gridflag_6;
int gewaldflag,gewaldflag_6;
int minorder,overlap_allowed;
int adjust_cutoff_flag;
int suffix_flag; // suffix compatibility flag
double scale;
double **gcons,**dgcons; // accumulated per-atom energy/virial
int evflag,evflag_atom;
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int maxeatom,maxvatom;
void pair_check();
void ev_setup(int, int);
double estimate_table_accuracy(double, double);
};
}
#endif
/* ERROR/WARNING messages:
E: KSpace solver requires a pair style
No pair style is defined.
E: KSpace style is incompatible with Pair style
Setting a kspace style requires that a pair style with a long-range
Coulombic or dispersion component be used.
W: For better accuracy use 'pair_modify table 0'
The user-specified force accuracy cannot be achieved unless the table
feature is disabled by using 'pair_modify table 0'.
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: Bad kspace_modify slab parameter
Kspace_modify value for the slab/volume keyword must be >= 2.0.
W: Kspace_modify slab param < 2.0 may cause unphysical behavior
The kspace_modify slab parameter should be larger to insure periodic
grids padded with empty space do not overlap.
*/