Merge branch 'master' of https://github.com/lammps/lammps into kk_hash
This commit is contained in:
Stan Moore
@ -25,6 +25,11 @@
|
||||
#include <string>
|
||||
#include <zstd.h>
|
||||
|
||||
#if ZSTD_VERSION_NUMBER < 10400
|
||||
#error must have at least zstd version 1.4 to compile with -DLAMMPS_ZSTD
|
||||
#endif
|
||||
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ZstdFileWriter : public FileWriter {
|
||||
|
||||
@ -31,7 +31,7 @@ action () {
|
||||
|
||||
if (test $1 = 1) then
|
||||
if (test ! -e ../pppm.cpp) then
|
||||
echo "Must install KSPACE package with DIELECTRIC"
|
||||
echo "Must install KSPACE package with DIELECTRIC package"
|
||||
exit 1
|
||||
fi
|
||||
fi
|
||||
|
||||
@ -431,7 +431,7 @@ void FixPolarizeBEMGMRES::gmres_solve(double *x, double *r)
|
||||
|
||||
// fill up h with zero
|
||||
|
||||
memset(h, 0, (mr + 1) * mr * sizeof(double));
|
||||
memset(h, 0, (size_t)(mr + 1) * mr * sizeof(double));
|
||||
|
||||
// the inner loop k = 1..(n-1)
|
||||
// build up the k-th Krylov space,
|
||||
@ -756,11 +756,11 @@ double FixPolarizeBEMGMRES::memory_usage()
|
||||
bytes += atom->nmax * sizeof(double); // q_backup
|
||||
bytes += mr * sizeof(double); // c
|
||||
bytes += (mr + 1) * sizeof(double); // g
|
||||
bytes += (mr + 1) * mr * sizeof(double); // h
|
||||
bytes += (double) (mr + 1) * mr * sizeof(double); // h
|
||||
bytes += mat_dim * sizeof(double); // r
|
||||
bytes += mr * (mr + 1) * sizeof(double); // s
|
||||
bytes += (double) mr * (mr + 1) * sizeof(double); // s
|
||||
bytes += mat_dim * sizeof(double); // v
|
||||
bytes += (mr + 1) * mr * sizeof(double); // y
|
||||
bytes += (double) (mr + 1) * mr * sizeof(double); // y
|
||||
return bytes;
|
||||
}
|
||||
|
||||
|
||||
@ -37,6 +37,7 @@
|
||||
#include "kspace.h"
|
||||
#include "math_const.h"
|
||||
#include "math_extra.h"
|
||||
#include "math_special.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "msm_dielectric.h"
|
||||
@ -60,6 +61,7 @@ using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathExtra;
|
||||
using namespace MathConst;
|
||||
using namespace MathSpecial;
|
||||
|
||||
enum { REAL2SCALED = 0, SCALED2REAL = 1 };
|
||||
|
||||
@ -589,21 +591,21 @@ void FixPolarizeFunctional::set_arrays(int i)
|
||||
double FixPolarizeFunctional::memory_usage()
|
||||
{
|
||||
double bytes = 0;
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // inverse_matrix
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // Rww
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G1ww
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // ndotGww
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G2ww
|
||||
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G3ww
|
||||
bytes += num_induced_charges * sizeof(double); // qiRqwVector
|
||||
bytes += num_induced_charges * sizeof(double); // sum2G2wq
|
||||
bytes += num_induced_charges * sizeof(double); // sum1G2qw
|
||||
bytes += num_induced_charges * sizeof(double); // sum1G1qw_epsilon
|
||||
bytes += num_induced_charges * sizeof(double); // sum2ndotGwq_epsilon
|
||||
bytes += num_ions * num_induced_charges * sizeof(double); // G1qw_real
|
||||
bytes += nmax * sizeof(int); // induced_charge_idx
|
||||
bytes += nmax * sizeof(int); // ion_idx
|
||||
bytes += num_induced_charges * sizeof(double); // induced_charges
|
||||
bytes += square(num_induced_charges) * sizeof(double); // inverse_matrix
|
||||
bytes += square(num_induced_charges) * sizeof(double); // Rww
|
||||
bytes += square(num_induced_charges) * sizeof(double); // G1ww
|
||||
bytes += square(num_induced_charges) * sizeof(double); // ndotGww
|
||||
bytes += square(num_induced_charges) * sizeof(double); // G2ww
|
||||
bytes += square(num_induced_charges) * sizeof(double); // G3ww
|
||||
bytes += num_induced_charges * sizeof(double); // qiRqwVector
|
||||
bytes += num_induced_charges * sizeof(double); // sum2G2wq
|
||||
bytes += num_induced_charges * sizeof(double); // sum1G2qw
|
||||
bytes += num_induced_charges * sizeof(double); // sum1G1qw_epsilon
|
||||
bytes += num_induced_charges * sizeof(double); // sum2ndotGwq_epsilon
|
||||
bytes += (double)num_ions * num_induced_charges * sizeof(double); // G1qw_real
|
||||
bytes += nmax * sizeof(int); // induced_charge_idx
|
||||
bytes += nmax * sizeof(int); // ion_idx
|
||||
bytes += num_induced_charges * sizeof(double); // induced_charges
|
||||
return bytes;
|
||||
}
|
||||
|
||||
@ -1115,9 +1117,6 @@ void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec, double dx,
|
||||
void FixPolarizeFunctional::calculate_matrix_multiply_vector(double **matrix, double *in_vec,
|
||||
double *out_vec, int M)
|
||||
{
|
||||
#if defined(OPENMP)
|
||||
#pragma parallel omp for
|
||||
#endif
|
||||
for (int k = 0; k < M; ++k) {
|
||||
double temp = 0.0;
|
||||
for (int l = 0; l < M; ++l) { temp += matrix[k][l] * in_vec[l]; }
|
||||
|
||||
@ -53,9 +53,9 @@ using namespace MathExtra;
|
||||
PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp) : PairLJLongCoulLong(lmp)
|
||||
{
|
||||
respa_enable = 0;
|
||||
cut_respa = NULL;
|
||||
efield = NULL;
|
||||
epot = NULL;
|
||||
cut_respa = nullptr;
|
||||
efield = nullptr;
|
||||
epot = nullptr;
|
||||
nmax = 0;
|
||||
}
|
||||
|
||||
|
||||
@ -469,21 +469,19 @@ void PPPMDielectric::slabcorr()
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute qsum,qsqsum,q2 and give error/warning if not charge neutral
|
||||
called initially, when particle count changes, when charges are changed
|
||||
compute qsum,qsqsum,q2 and ignore error/warning if not charge neutral
|
||||
called whenever charges are changed
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDielectric::qsum_qsq()
|
||||
{
|
||||
const double * const q = atom->q;
|
||||
const double * const eps = atom->epsilon;
|
||||
const int nlocal = atom->nlocal;
|
||||
double qsum_local(0.0), qsqsum_local(0.0);
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
double qtmp = eps[i]*q[i];
|
||||
qsum_local += qtmp;
|
||||
qsqsum_local += qtmp*qtmp;
|
||||
qsum_local += q[i];
|
||||
qsqsum_local += q[i]*q[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
868
src/DIELECTRIC/pppm_disp_dielectric.cpp
Normal file
868
src/DIELECTRIC/pppm_disp_dielectric.cpp
Normal file
@ -0,0 +1,868 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ 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: Trung Nguyen (Northwestern)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pppm_disp_dielectric.h"
|
||||
|
||||
#include "angle.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_dielectric.h"
|
||||
#include "bond.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fft3d_wrap.h"
|
||||
#include "force.h"
|
||||
#include "gridcomm.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "remap_wrap.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
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,REVERSE_RHO_GEOM,REVERSE_RHO_ARITH,REVERSE_RHO_NONE};
|
||||
enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM,
|
||||
FORWARD_IK_GEOM,FORWARD_AD_GEOM,
|
||||
FORWARD_IK_PERATOM_GEOM,FORWARD_AD_PERATOM_GEOM,
|
||||
FORWARD_IK_ARITH,FORWARD_AD_ARITH,
|
||||
FORWARD_IK_PERATOM_ARITH,FORWARD_AD_PERATOM_ARITH,
|
||||
FORWARD_IK_NONE,FORWARD_AD_NONE,FORWARD_IK_PERATOM_NONE,
|
||||
FORWARD_AD_PERATOM_NONE};
|
||||
|
||||
#ifdef FFT_SINGLE
|
||||
#define ZEROF 0.0f
|
||||
#define ONEF 1.0f
|
||||
#else
|
||||
#define ZEROF 0.0
|
||||
#define ONEF 1.0
|
||||
#endif
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PPPMDispDielectric::PPPMDispDielectric(LAMMPS *lmp) : PPPMDisp(lmp)
|
||||
{
|
||||
dipoleflag = 0; // turned off for now, until dipole works
|
||||
group_group_enable = 0;
|
||||
|
||||
mu_flag = 0;
|
||||
|
||||
efield = nullptr;
|
||||
phi = nullptr;
|
||||
potflag = 0;
|
||||
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"pppm/dielectric requires atom style dielectric");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PPPMDispDielectric::~PPPMDispDielectric()
|
||||
{
|
||||
memory->destroy(efield);
|
||||
memory->destroy(phi);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute the PPPM long-range force, energy, virial
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
|
||||
int i;
|
||||
|
||||
// set energy/virial flags
|
||||
// invoke allocate_peratom() if needed for first time
|
||||
|
||||
ev_init(eflag,vflag);
|
||||
|
||||
if (evflag_atom && !peratom_allocate_flag) allocate_peratom();
|
||||
|
||||
// convert atoms from box to lamda coords
|
||||
|
||||
if (triclinic == 0) boxlo = domain->boxlo;
|
||||
else {
|
||||
boxlo = domain->boxlo_lamda;
|
||||
domain->x2lamda(atom->nlocal);
|
||||
}
|
||||
// extend size of per-atom arrays if necessary
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
|
||||
if (function[0]) {
|
||||
memory->destroy(part2grid);
|
||||
memory->destroy(efield);
|
||||
memory->destroy(phi);
|
||||
}
|
||||
if (function[1] + function[2] + function[3]) memory->destroy(part2grid_6);
|
||||
nmax = atom->nmax;
|
||||
if (function[0]) {
|
||||
memory->create(part2grid,nmax,3,"pppm/disp:part2grid");
|
||||
memory->create(efield,nmax,3,"pppm/disp:efield");
|
||||
memory->create(phi,nmax,"pppm/disp:phi");
|
||||
}
|
||||
if (function[1] + function[2] + function[3])
|
||||
memory->create(part2grid_6,nmax,3,"pppm/disp:part2grid_6");
|
||||
}
|
||||
|
||||
energy = 0.0;
|
||||
energy_1 = 0.0;
|
||||
energy_6 = 0.0;
|
||||
if (vflag) for (i = 0; i < 6; i++) virial_6[i] = virial_1[i] = 0.0;
|
||||
|
||||
// find grid points for all my particles
|
||||
// distribute partcles' charges/dispersion coefficients on the grid
|
||||
// communication between processors and remapping two fft
|
||||
// Solution of poissons equation in k-space and backtransformation
|
||||
// communication between processors
|
||||
// calculation of forces
|
||||
|
||||
if (function[0]) {
|
||||
|
||||
// perform calculations for coulomb interactions only
|
||||
|
||||
particle_map_c(delxinv,delyinv,delzinv,shift,part2grid,nupper,nlower,
|
||||
nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out);
|
||||
|
||||
make_rho_c();
|
||||
|
||||
gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO,
|
||||
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
brick2fft(nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
|
||||
density_brick,density_fft,work1,remap);
|
||||
|
||||
if (differentiation_flag == 1) {
|
||||
poisson_ad(work1,work2,density_fft,fft1,fft2,
|
||||
nx_pppm,ny_pppm,nz_pppm,nfft,
|
||||
nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
|
||||
nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
|
||||
energy_1,greensfn,
|
||||
virial_1,vg,vg2,
|
||||
u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
|
||||
|
||||
gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD,
|
||||
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_c_ad();
|
||||
|
||||
if (vflag_atom)
|
||||
gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM,
|
||||
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
} else {
|
||||
poisson_ik(work1,work2,density_fft,fft1,fft2,
|
||||
nx_pppm,ny_pppm,nz_pppm,nfft,
|
||||
nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
|
||||
nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
|
||||
energy_1,greensfn,
|
||||
fkx,fky,fkz,fkx2,fky2,fkz2,
|
||||
vdx_brick,vdy_brick,vdz_brick,virial_1,vg,vg2,
|
||||
u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
|
||||
|
||||
gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK,
|
||||
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_c_ik();
|
||||
|
||||
if (evflag_atom)
|
||||
gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM,
|
||||
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
|
||||
}
|
||||
|
||||
if (evflag_atom) fieldforce_c_peratom();
|
||||
}
|
||||
|
||||
if (function[1]) {
|
||||
|
||||
// perform calculations for geometric mixing
|
||||
|
||||
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
|
||||
nupper_6,nlower_6,
|
||||
nxlo_out_6,nylo_out_6,nzlo_out_6,
|
||||
nxhi_out_6,nyhi_out_6,nzhi_out_6);
|
||||
|
||||
make_rho_g();
|
||||
|
||||
gc6->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_GEOM,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
brick2fft(nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
|
||||
density_brick_g,density_fft_g,work1_6,remap_6);
|
||||
|
||||
if (differentiation_flag == 1) {
|
||||
poisson_ad(work1_6,work2_6,density_fft_g,fft1_6,fft2_6,
|
||||
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
|
||||
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
|
||||
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
|
||||
energy_6,greensfn_6,
|
||||
virial_6,vg_6,vg2_6,
|
||||
u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g,
|
||||
v3_brick_g,v4_brick_g,v5_brick_g);
|
||||
|
||||
gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_GEOM,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_g_ad();
|
||||
|
||||
if (vflag_atom)
|
||||
gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_GEOM,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
} else {
|
||||
poisson_ik(work1_6,work2_6,density_fft_g,fft1_6,fft2_6,
|
||||
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
|
||||
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
|
||||
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
|
||||
energy_6,greensfn_6,
|
||||
fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6,
|
||||
vdx_brick_g,vdy_brick_g,vdz_brick_g,virial_6,vg_6,vg2_6,
|
||||
u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g,
|
||||
v3_brick_g,v4_brick_g,v5_brick_g);
|
||||
|
||||
gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_GEOM,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_g_ik();
|
||||
|
||||
if (evflag_atom)
|
||||
gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_GEOM,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
}
|
||||
|
||||
if (evflag_atom) fieldforce_g_peratom();
|
||||
}
|
||||
|
||||
if (function[2]) {
|
||||
|
||||
// perform calculations for arithmetic mixing
|
||||
|
||||
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
|
||||
nupper_6,nlower_6,
|
||||
nxlo_out_6,nylo_out_6,nzlo_out_6,
|
||||
nxhi_out_6,nyhi_out_6,nzhi_out_6);
|
||||
|
||||
make_rho_a();
|
||||
|
||||
gc6->reverse_comm_kspace(this,7,sizeof(FFT_SCALAR),REVERSE_RHO_ARITH,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
brick2fft_a();
|
||||
|
||||
if (differentiation_flag == 1) {
|
||||
poisson_ad(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6,
|
||||
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
|
||||
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
|
||||
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
|
||||
energy_6,greensfn_6,
|
||||
virial_6,vg_6,vg2_6,
|
||||
u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3,
|
||||
v3_brick_a3,v4_brick_a3,v5_brick_a3);
|
||||
poisson_2s_ad(density_fft_a0,density_fft_a6,
|
||||
u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0,
|
||||
v3_brick_a0,v4_brick_a0,v5_brick_a0,
|
||||
u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6,
|
||||
v3_brick_a6,v4_brick_a6,v5_brick_a6);
|
||||
poisson_2s_ad(density_fft_a1,density_fft_a5,
|
||||
u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1,
|
||||
v3_brick_a1,v4_brick_a1,v5_brick_a1,
|
||||
u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5,
|
||||
v3_brick_a5,v4_brick_a5,v5_brick_a5);
|
||||
poisson_2s_ad(density_fft_a2,density_fft_a4,
|
||||
u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2,
|
||||
v3_brick_a2,v4_brick_a2,v5_brick_a2,
|
||||
u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4,
|
||||
v3_brick_a4,v4_brick_a4,v5_brick_a4);
|
||||
|
||||
gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_ARITH,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_a_ad();
|
||||
|
||||
if (evflag_atom)
|
||||
gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_ARITH,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
} else {
|
||||
poisson_ik(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6,
|
||||
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
|
||||
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
|
||||
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
|
||||
energy_6,greensfn_6,
|
||||
fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6,
|
||||
vdx_brick_a3,vdy_brick_a3,vdz_brick_a3,virial_6,vg_6,vg2_6,
|
||||
u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3,
|
||||
v3_brick_a3,v4_brick_a3,v5_brick_a3);
|
||||
poisson_2s_ik(density_fft_a0,density_fft_a6,
|
||||
vdx_brick_a0,vdy_brick_a0,vdz_brick_a0,
|
||||
vdx_brick_a6,vdy_brick_a6,vdz_brick_a6,
|
||||
u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0,
|
||||
v3_brick_a0,v4_brick_a0,v5_brick_a0,
|
||||
u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6,
|
||||
v3_brick_a6,v4_brick_a6,v5_brick_a6);
|
||||
poisson_2s_ik(density_fft_a1,density_fft_a5,
|
||||
vdx_brick_a1,vdy_brick_a1,vdz_brick_a1,
|
||||
vdx_brick_a5,vdy_brick_a5,vdz_brick_a5,
|
||||
u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1,
|
||||
v3_brick_a1,v4_brick_a1,v5_brick_a1,
|
||||
u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5,
|
||||
v3_brick_a5,v4_brick_a5,v5_brick_a5);
|
||||
poisson_2s_ik(density_fft_a2,density_fft_a4,
|
||||
vdx_brick_a2,vdy_brick_a2,vdz_brick_a2,
|
||||
vdx_brick_a4,vdy_brick_a4,vdz_brick_a4,
|
||||
u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2,
|
||||
v3_brick_a2,v4_brick_a2,v5_brick_a2,
|
||||
u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4,
|
||||
v3_brick_a4,v4_brick_a4,v5_brick_a4);
|
||||
|
||||
gc6->forward_comm_kspace(this,21,sizeof(FFT_SCALAR),FORWARD_IK_ARITH,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_a_ik();
|
||||
|
||||
if (evflag_atom)
|
||||
gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_ARITH,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
}
|
||||
|
||||
if (evflag_atom) fieldforce_a_peratom();
|
||||
}
|
||||
|
||||
if (function[3]) {
|
||||
|
||||
// perform calculations if no mixing rule applies
|
||||
|
||||
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
|
||||
nupper_6,nlower_6,
|
||||
nxlo_out_6,nylo_out_6,nzlo_out_6,
|
||||
nxhi_out_6,nyhi_out_6,nzhi_out_6);
|
||||
|
||||
make_rho_none();
|
||||
|
||||
gc6->reverse_comm_kspace(this,nsplit_alloc,sizeof(FFT_SCALAR),REVERSE_RHO_NONE,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
brick2fft_none();
|
||||
|
||||
if (differentiation_flag == 1) {
|
||||
int n = 0;
|
||||
for (int k = 0; k < nsplit_alloc/2; k++) {
|
||||
poisson_none_ad(n,n+1,density_fft_none[n],density_fft_none[n+1],
|
||||
u_brick_none[n],u_brick_none[n+1],
|
||||
v0_brick_none,v1_brick_none,v2_brick_none,
|
||||
v3_brick_none,v4_brick_none,v5_brick_none);
|
||||
n += 2;
|
||||
}
|
||||
|
||||
gc6->forward_comm_kspace(this,1*nsplit_alloc,sizeof(FFT_SCALAR),
|
||||
FORWARD_AD_NONE,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_none_ad();
|
||||
|
||||
if (vflag_atom)
|
||||
gc6->forward_comm_kspace(this,6*nsplit_alloc,sizeof(FFT_SCALAR),
|
||||
FORWARD_AD_PERATOM_NONE,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
} else {
|
||||
int n = 0;
|
||||
for (int k = 0; k < nsplit_alloc/2; k++) {
|
||||
poisson_none_ik(n,n+1,density_fft_none[n],density_fft_none[n+1],
|
||||
vdx_brick_none[n],vdy_brick_none[n],vdz_brick_none[n],
|
||||
vdx_brick_none[n+1],vdy_brick_none[n+1],vdz_brick_none[n+1],
|
||||
u_brick_none,v0_brick_none,v1_brick_none,v2_brick_none,
|
||||
v3_brick_none,v4_brick_none,v5_brick_none);
|
||||
n += 2;
|
||||
}
|
||||
|
||||
gc6->forward_comm_kspace(this,3*nsplit_alloc,sizeof(FFT_SCALAR),
|
||||
FORWARD_IK_NONE,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
|
||||
fieldforce_none_ik();
|
||||
|
||||
if (evflag_atom)
|
||||
gc6->forward_comm_kspace(this,7*nsplit_alloc,sizeof(FFT_SCALAR),
|
||||
FORWARD_IK_PERATOM_NONE,
|
||||
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
|
||||
}
|
||||
|
||||
if (evflag_atom) fieldforce_none_peratom();
|
||||
}
|
||||
|
||||
// update qsum and qsqsum, if atom count has changed and energy needed
|
||||
|
||||
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
|
||||
qsum_qsq();
|
||||
natoms_original = atom->natoms;
|
||||
}
|
||||
|
||||
// sum energy across procs and add in volume-dependent term
|
||||
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
|
||||
if (eflag_global) {
|
||||
double energy_all;
|
||||
MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
energy_1 = energy_all;
|
||||
MPI_Allreduce(&energy_6,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
energy_6 = energy_all;
|
||||
|
||||
energy_1 *= 0.5*volume;
|
||||
energy_6 *= 0.5*volume;
|
||||
|
||||
energy_1 -= g_ewald*qsqsum/MY_PIS +
|
||||
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
|
||||
energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +
|
||||
1.0/12.0*pow(g_ewald_6,6)*csum;
|
||||
energy_1 *= qscale;
|
||||
}
|
||||
|
||||
// sum virial across procs
|
||||
|
||||
if (vflag_global) {
|
||||
double virial_all[6];
|
||||
MPI_Allreduce(virial_1,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
|
||||
MPI_Allreduce(virial_6,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
for (i = 0; i < 6; i++) virial[i] += 0.5*volume*virial_all[i];
|
||||
if (function[1]+function[2]+function[3]) {
|
||||
double a = MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij;
|
||||
virial[0] -= a;
|
||||
virial[1] -= a;
|
||||
virial[2] -= a;
|
||||
}
|
||||
}
|
||||
|
||||
if (eflag_atom) {
|
||||
if (function[0]) {
|
||||
double *q = atom->q;
|
||||
// coulomb self energy correction
|
||||
for (i = 0; i < atom->nlocal; i++) {
|
||||
eatom[i] -= qscale*g_ewald*q[i]*q[i]/MY_PIS +
|
||||
qscale*MY_PI2*q[i]*qsum / (g_ewald*g_ewald*volume);
|
||||
}
|
||||
}
|
||||
if (function[1] + function[2] + function[3]) {
|
||||
int tmp;
|
||||
for (i = 0; i < atom->nlocal; i++) {
|
||||
tmp = atom->type[i];
|
||||
eatom[i] += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp] +
|
||||
1.0/12.0*pow(g_ewald_6,6)*cii[tmp];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_atom) {
|
||||
if (function[1] + function[2] + function[3]) {
|
||||
int tmp;
|
||||
// dispersion self virial correction
|
||||
for (i = 0; i < atom->nlocal; i++) {
|
||||
tmp = atom->type[i];
|
||||
for (int n = 0; n < 3; n++)
|
||||
vatom[i][n] -= MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 2d slab correction
|
||||
|
||||
if (slabflag) slabcorr(eflag);
|
||||
if (function[0]) energy += energy_1;
|
||||
if (function[1] + function[2] + function[3]) energy += energy_6;
|
||||
|
||||
// convert atoms back from lamda to box coords
|
||||
|
||||
if (triclinic) domain->lamda2x(atom->nlocal);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get electric field & force on my particles
|
||||
for ik scheme
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::fieldforce_c_ik()
|
||||
{
|
||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
||||
FFT_SCALAR ekx,eky,ekz,u;
|
||||
|
||||
// loop over my charges, interpolate electric field from nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
// ek = 3 components of E-field on particle
|
||||
|
||||
double *q = atom->q;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *eps = atom->epsilon;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
|
||||
|
||||
u = ekx = eky = ekz = ZEROF;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
z0 = rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
y0 = z0*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
x0 = y0*rho1d[0][l];
|
||||
if (potflag) u += x0*u_brick[mz][my][mx];
|
||||
ekx -= x0*vdx_brick[mz][my][mx];
|
||||
eky -= x0*vdy_brick[mz][my][mx];
|
||||
ekz -= x0*vdz_brick[mz][my][mx];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// electrostatic potential
|
||||
|
||||
if (potflag) phi[i] = u;
|
||||
|
||||
// convert E-field to force
|
||||
|
||||
const double efactor = scale * eps[i];
|
||||
efield[i][0] = efactor*ekx;
|
||||
efield[i][1] = efactor*eky;
|
||||
efield[i][2] = efactor*ekz;
|
||||
|
||||
// convert E-field to force
|
||||
|
||||
const double qfactor = force->qqrd2e * scale * q[i];
|
||||
f[i][0] += qfactor*ekx;
|
||||
f[i][1] += qfactor*eky;
|
||||
if (slabflag != 2) f[i][2] += qfactor*ekz;
|
||||
}
|
||||
}
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get electric field & force on my particles
|
||||
for ad scheme
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::fieldforce_c_ad()
|
||||
{
|
||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz;
|
||||
FFT_SCALAR ekx,eky,ekz,u;
|
||||
double s1,s2,s3;
|
||||
double sf = 0.0;
|
||||
|
||||
double *prd;
|
||||
|
||||
if (triclinic == 0) prd = domain->prd;
|
||||
else prd = domain->prd_lamda;
|
||||
|
||||
double xprd = prd[0];
|
||||
double yprd = prd[1];
|
||||
double zprd = prd[2];
|
||||
double zprd_slab = zprd*slab_volfactor;
|
||||
|
||||
double hx_inv = nx_pppm/xprd;
|
||||
double hy_inv = ny_pppm/yprd;
|
||||
double hz_inv = nz_pppm/zprd_slab;
|
||||
|
||||
// loop over my charges, interpolate electric field from nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
// ek = 3 components of E-field on particle
|
||||
|
||||
double *q = atom->q;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *eps = atom->epsilon;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
|
||||
compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d);
|
||||
|
||||
u = ekx = eky = ekz = ZEROF;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
u += rho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
|
||||
ekx += drho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
|
||||
eky += rho1d[0][l]*drho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
|
||||
ekz += rho1d[0][l]*rho1d[1][m]*drho1d[2][n]*u_brick[mz][my][mx];
|
||||
}
|
||||
}
|
||||
}
|
||||
ekx *= hx_inv;
|
||||
eky *= hy_inv;
|
||||
ekz *= hz_inv;
|
||||
|
||||
// electrical potential
|
||||
|
||||
if (potflag) phi[i] = u;
|
||||
|
||||
// convert E-field to force and substract self forces
|
||||
const double qfactor = qqrd2e * scale;
|
||||
double qtmp = eps[i]*q[i];
|
||||
|
||||
s1 = x[i][0]*hx_inv;
|
||||
s2 = x[i][1]*hy_inv;
|
||||
s3 = x[i][2]*hz_inv;
|
||||
sf = sf_coeff[0]*sin(2*MY_PI*s1);
|
||||
sf += sf_coeff[1]*sin(4*MY_PI*s1);
|
||||
sf *= 2*q[i]*q[i];
|
||||
f[i][0] += qfactor*(ekx*q[i] - sf);
|
||||
|
||||
sf = sf_coeff[2]*sin(2*MY_PI*s2);
|
||||
sf += sf_coeff[3]*sin(4*MY_PI*s2);
|
||||
sf *= 2*q[i]*q[i];
|
||||
f[i][1] += qfactor*(eky*q[i] - sf);
|
||||
|
||||
sf = sf_coeff[4]*sin(2*MY_PI*s3);
|
||||
sf += sf_coeff[5]*sin(4*MY_PI*s3);
|
||||
sf *= 2*q[i]*q[i];
|
||||
if (slabflag != 2) f[i][2] += qfactor*(ekz*q[i] - sf);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get electric field & force on my particles
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::fieldforce_c_peratom()
|
||||
{
|
||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
||||
FFT_SCALAR u_pa,v0,v1,v2,v3,v4,v5;
|
||||
|
||||
// loop over my charges, interpolate electric field from nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
// ek = 3 components of E-field on particle
|
||||
|
||||
double *q = atom->q;
|
||||
double **x = atom->x;
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
|
||||
|
||||
u_pa = v0 = v1 = v2 = v3 = v4 = v5 = ZEROF;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
z0 = rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
y0 = z0*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
x0 = y0*rho1d[0][l];
|
||||
if (eflag_atom) u_pa += x0*u_brick[mz][my][mx];
|
||||
if (vflag_atom) {
|
||||
v0 += x0*v0_brick[mz][my][mx];
|
||||
v1 += x0*v1_brick[mz][my][mx];
|
||||
v2 += x0*v2_brick[mz][my][mx];
|
||||
v3 += x0*v3_brick[mz][my][mx];
|
||||
v4 += x0*v4_brick[mz][my][mx];
|
||||
v5 += x0*v5_brick[mz][my][mx];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// electrostatic potential
|
||||
|
||||
phi[i] = u_pa;
|
||||
|
||||
// convert E-field to force
|
||||
|
||||
const double qfactor = 0.5*force->qqrd2e * scale * q[i];
|
||||
|
||||
if (eflag_atom) eatom[i] += u_pa*qfactor;
|
||||
if (vflag_atom) {
|
||||
vatom[i][0] += v0*qfactor;
|
||||
vatom[i][1] += v1*qfactor;
|
||||
vatom[i][2] += v2*qfactor;
|
||||
vatom[i][3] += v3*qfactor;
|
||||
vatom[i][4] += v4*qfactor;
|
||||
vatom[i][5] += v5*qfactor;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Slab-geometry correction term to dampen inter-slab interactions between
|
||||
periodically repeating slabs. Yields good approximation to 2D Ewald if
|
||||
adequate empty space is left between repeating slabs (J. Chem. Phys.
|
||||
111, 3155). Slabs defined here to be parallel to the xy plane. Also
|
||||
extended to non-neutral systems (J. Chem. Phys. 131, 094107).
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::slabcorr(int eflag)
|
||||
{
|
||||
// compute local contribution to global dipole moment
|
||||
|
||||
double *q = atom->q;
|
||||
double **x = atom->x;
|
||||
double *eps = atom->epsilon;
|
||||
double zprd = domain->zprd;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double dipole = 0.0;
|
||||
for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
|
||||
|
||||
if (mu_flag) {
|
||||
double **mu = atom->mu;
|
||||
for (int i = 0; i < nlocal; i++) dipole += mu[i][2];
|
||||
}
|
||||
|
||||
// sum local contributions to get global dipole moment
|
||||
|
||||
double dipole_all;
|
||||
MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// need to make non-neutral systems and/or
|
||||
// per-atom energy translationally invariant
|
||||
|
||||
double dipole_r2 = 0.0;
|
||||
if (eflag_atom || fabs(qsum) > SMALL) {
|
||||
|
||||
if (mu_flag)
|
||||
error->all(FLERR,"Cannot (yet) use kspace slab correction with "
|
||||
"long-range dipoles and non-neutral systems or per-atom energy");
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
dipole_r2 += q[i]*x[i][2]*x[i][2];
|
||||
|
||||
// sum local contributions
|
||||
|
||||
double tmp;
|
||||
MPI_Allreduce(&dipole_r2,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
dipole_r2 = tmp;
|
||||
}
|
||||
|
||||
// compute corrections
|
||||
|
||||
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
|
||||
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
|
||||
const double qscale = qqrd2e * scale;
|
||||
|
||||
if (eflag_global) energy += qscale * e_slabcorr;
|
||||
|
||||
// per-atom energy
|
||||
|
||||
if (eflag_atom) {
|
||||
double efact = qscale * MY_2PI/volume;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
eatom[i] += efact * eps[i]*q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 +
|
||||
qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0);
|
||||
}
|
||||
|
||||
// add on force corrections
|
||||
|
||||
double ffact = qscale * (-4.0*MY_PI/volume);
|
||||
double **f = atom->f;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
f[i][2] += ffact * eps[i]*q[i]*(dipole_all - qsum*x[i][2]);
|
||||
efield[i][2] += ffact * eps[i]*(dipole_all - qsum*x[i][2]);
|
||||
}
|
||||
|
||||
// add on torque corrections
|
||||
|
||||
if (mu_flag && atom->torque) {
|
||||
double **mu = atom->mu;
|
||||
double **torque = atom->torque;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
torque[i][0] += ffact * dipole_all * mu[i][1];
|
||||
torque[i][1] += -ffact * dipole_all * mu[i][0];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PPPMDispDielectric::memory_usage()
|
||||
{
|
||||
double bytes = PPPMDisp::memory_usage();
|
||||
bytes += nmax*3 * sizeof(double);
|
||||
bytes += nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute qsum,qsqsum,q2 and give error/warning if not charge neutral
|
||||
called initially, when particle count changes, when charges are changed
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDispDielectric::qsum_qsq()
|
||||
{
|
||||
const double * const q = atom->q;
|
||||
const int nlocal = atom->nlocal;
|
||||
double qsum_local(0.0), qsqsum_local(0.0);
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
qsum_local += q[i];
|
||||
qsqsum_local += q[i]*q[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
q2 = qsqsum * force->qqrd2e;
|
||||
}
|
||||
62
src/DIELECTRIC/pppm_disp_dielectric.h
Normal file
62
src/DIELECTRIC/pppm_disp_dielectric.h
Normal file
@ -0,0 +1,62 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ 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
|
||||
// clang-format off
|
||||
KSpaceStyle(pppm/disp/dielectric,PPPMDispDielectric);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PPPM_DISP_DIELECTRIC_H
|
||||
#define LMP_PPPM_DISP_DIELECTRIC_H
|
||||
|
||||
#include "pppm_disp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PPPMDispDielectric : public PPPMDisp {
|
||||
public:
|
||||
PPPMDispDielectric(class LAMMPS *);
|
||||
virtual ~PPPMDispDielectric();
|
||||
virtual double memory_usage();
|
||||
virtual void compute(int, int);
|
||||
void qsum_qsq();
|
||||
void slabcorr(int);
|
||||
|
||||
double **efield;
|
||||
double *phi;
|
||||
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
|
||||
|
||||
protected:
|
||||
virtual void fieldforce_c_ik();
|
||||
virtual void fieldforce_c_ad();
|
||||
virtual void fieldforce_c_peratom();
|
||||
|
||||
class AtomVecDielectric *avec;
|
||||
int mu_flag;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#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.
|
||||
|
||||
*/
|
||||
@ -78,6 +78,7 @@ if (test $1 = "GRANULAR") then
|
||||
fi
|
||||
|
||||
if (test $1 = "KSPACE") then
|
||||
depend CG-SDK
|
||||
depend CORESHELL
|
||||
depend DIELECTRIC
|
||||
depend GPU
|
||||
|
||||
@ -610,7 +610,9 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
|
||||
|
||||
if (test) return;
|
||||
#else
|
||||
dev.team_barrier();
|
||||
int not_done = (i >= 0 && i <= nlocal);
|
||||
dev.team_reduce(Kokkos::Max<int>(not_done));
|
||||
if(not_done == 0) return;
|
||||
#endif
|
||||
|
||||
if (i >= 0 && i < nlocal) {
|
||||
@ -1053,13 +1055,14 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
|
||||
other_x[MY_II + 4 * atoms_per_bin] = radi;
|
||||
}
|
||||
other_id[MY_II] = i;
|
||||
// FIXME_SYCL
|
||||
#ifndef KOKKOS_ENABLE_SYCL
|
||||
int test = (__syncthreads_count(i >= 0 && i <= nlocal) == 0);
|
||||
|
||||
if (test) return;
|
||||
#else
|
||||
dev.team_barrier();
|
||||
int not_done = (i >= 0 && i <= nlocal);
|
||||
dev.team_reduce(Kokkos::Max<int>(not_done));
|
||||
if(not_done == 0) return;
|
||||
#endif
|
||||
|
||||
if (i >= 0 && i < nlocal) {
|
||||
|
||||
@ -705,7 +705,7 @@ void MLIAP_SO3::get_sbes_array(int nlocal, int *numneighs, double **rij, int lma
|
||||
bigint ipair = 0;
|
||||
bigint gindex;
|
||||
int findex = m_Nmax * (m_lmax + 1);
|
||||
int mindex = m_lmax + 1;
|
||||
const bigint mindex = m_lmax + 1;
|
||||
|
||||
for (int ii = 0; ii < nlocal; ii++) {
|
||||
|
||||
@ -721,10 +721,10 @@ void MLIAP_SO3::get_sbes_array(int nlocal, int *numneighs, double **rij, int lma
|
||||
if (ri < SMALL) continue;
|
||||
|
||||
pfac2 = pfac1 * ri;
|
||||
|
||||
gindex = (ipair - 1) * findex;
|
||||
|
||||
for (i = 1; i < m_Nmax + 1; i++) {
|
||||
const bigint i1mindex = (bigint) (i - 1) * mindex;
|
||||
|
||||
x = cos((2 * i - 1) * pfac3);
|
||||
xi = pfac4 * (x + 1);
|
||||
@ -733,28 +733,27 @@ void MLIAP_SO3::get_sbes_array(int nlocal, int *numneighs, double **rij, int lma
|
||||
sa = sinh(rb) / rb;
|
||||
sb = (cosh(rb) - sa) / rb;
|
||||
|
||||
m_sbes_array[gindex + (i - 1) * mindex + 0] = sa;
|
||||
m_sbes_array[gindex + (i - 1) * mindex + 1] = sb;
|
||||
m_sbes_array[gindex + i1mindex + 0] = sa;
|
||||
m_sbes_array[gindex + i1mindex + 1] = sb;
|
||||
|
||||
for (j = 2; j < lmax + 1; j++)
|
||||
m_sbes_array[gindex + (i - 1) * mindex + j] =
|
||||
m_sbes_array[gindex + (i - 1) * mindex + j - 2] -
|
||||
(2 * j - 1) / rb * m_sbes_array[gindex + (i - 1) * mindex + j - 1];
|
||||
m_sbes_array[gindex + i1mindex + j] = m_sbes_array[gindex + i1mindex + j - 2] -
|
||||
(2 * j - 1) / rb * m_sbes_array[gindex + i1mindex + j - 1];
|
||||
|
||||
exts = m_sbes_array[gindex + (i - 1) * mindex + j - 2] -
|
||||
(2 * j - 1) / rb * m_sbes_array[gindex + (i - 1) * mindex + j - 1];
|
||||
exts = m_sbes_array[gindex + i1mindex + j - 2] -
|
||||
(2 * j - 1) / rb * m_sbes_array[gindex + i1mindex + j - 1];
|
||||
|
||||
m_sbes_darray[gindex + (i - 1) * mindex + 0] = sb;
|
||||
m_sbes_darray[gindex + i1mindex + 0] = sb;
|
||||
|
||||
for (j = 1; j < lmax; j++)
|
||||
m_sbes_darray[gindex + (i - 1) * mindex + j] = xi *
|
||||
(j * m_sbes_array[gindex + (i - 1) * mindex + j - 1] +
|
||||
(j + 1) * m_sbes_array[gindex + (i - 1) * mindex + j + 1]) /
|
||||
m_sbes_darray[gindex + i1mindex + j] = xi *
|
||||
(j * m_sbes_array[gindex + i1mindex + j - 1] +
|
||||
(j + 1) * m_sbes_array[gindex + i1mindex + j + 1]) /
|
||||
(2 * j + 1);
|
||||
|
||||
m_sbes_darray[gindex + (i - 1) * mindex + j] = xi *
|
||||
(j * m_sbes_array[gindex + (i - 1) * mindex + j - 1] + (j + 1) * exts) / (2 * j + 1);
|
||||
m_sbes_darray[gindex + (i - 1) * mindex + 0] = xi * sb;
|
||||
m_sbes_darray[gindex + i1mindex + j] =
|
||||
xi * (j * m_sbes_array[gindex + i1mindex + j - 1] + (j + 1) * exts) / (2 * j + 1);
|
||||
m_sbes_darray[gindex + i1mindex + 0] = xi * sb;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -797,14 +796,14 @@ void MLIAP_SO3::get_rip_array(int nlocal, int *numneighs, double **rij, int nmax
|
||||
integrald = 0.0;
|
||||
for (i = 0; i < m_Nmax; i++) {
|
||||
integral += m_g_array[(n - 1) * m_Nmax + i] *
|
||||
m_sbes_array[(ipair - 1) * m_Nmax * (m_lmax + 1) + i * (m_lmax + 1) + l];
|
||||
m_sbes_array[(ipair - 1) * m_Nmax * (m_lmax + 1) + (bigint) i * (m_lmax + 1) + l];
|
||||
integrald += m_g_array[(n - 1) * m_Nmax + i] *
|
||||
m_sbes_darray[(ipair - 1) * m_Nmax * (m_lmax + 1) + i * (m_lmax + 1) + l];
|
||||
m_sbes_darray[(ipair - 1) * m_Nmax * (m_lmax + 1) + (bigint) i * (m_lmax + 1) + l];
|
||||
}
|
||||
|
||||
m_rip_array[(ipair - 1) * m_nmax * (m_lmax + 1) + (n - 1) * (m_lmax + 1) + l] =
|
||||
m_rip_array[(ipair - 1) * m_nmax * (m_lmax + 1) + (bigint) (n - 1) * (m_lmax + 1) + l] =
|
||||
integral * expfac;
|
||||
m_rip_darray[(ipair - 1) * m_nmax * (m_lmax + 1) + (n - 1) * (m_lmax + 1) + l] =
|
||||
m_rip_darray[(ipair - 1) * m_nmax * (m_lmax + 1) + (bigint) (n - 1) * (m_lmax + 1) + l] =
|
||||
integrald * expfac;
|
||||
}
|
||||
}
|
||||
@ -904,7 +903,7 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem
|
||||
for (int n = 1; n < nmax + 1; n++) {
|
||||
int i = 0;
|
||||
for (int l = 0; l < lmax + 1; l++) {
|
||||
r_int = m_rip_array[gindex + (n - 1) * (m_lmax + 1) + l];
|
||||
r_int = m_rip_array[gindex + (bigint) (n - 1) * (m_lmax + 1) + l];
|
||||
|
||||
for (int m = -l; m < l + 1; m++) {
|
||||
|
||||
@ -1005,7 +1004,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
|
||||
|
||||
for (int ii = 0; ii < nlocal; ii++) {
|
||||
|
||||
totali = nmax * m_numYlms;
|
||||
totali = (bigint) nmax * m_numYlms;
|
||||
for (bigint ti = 0; ti < totali; ti++) {
|
||||
m_clisttot_r[ti] = 0.0;
|
||||
m_clisttot_i[ti] = 0.0;
|
||||
@ -1024,7 +1023,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
|
||||
r = sqrt(x * x + y * y + z * z);
|
||||
|
||||
if (r < SMALL) continue;
|
||||
totali = nmax * m_numYlms;
|
||||
totali = (bigint) nmax * m_numYlms;
|
||||
|
||||
for (bigint ti = 0; ti < totali; ti++) {
|
||||
m_clist_r[ti] = 0.0;
|
||||
@ -1044,7 +1043,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
|
||||
for (int n = 1; n < nmax + 1; n++) {
|
||||
int i = 0;
|
||||
for (int l = 0; l < lmax + 1; l++) {
|
||||
r_int = m_rip_array[gindex + (n - 1) * (m_lmax + 1) + l];
|
||||
r_int = m_rip_array[gindex + (bigint) (n - 1) * (m_lmax + 1) + l];
|
||||
|
||||
for (int m = -l; m < l + 1; m++) {
|
||||
|
||||
@ -1057,7 +1056,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
|
||||
}
|
||||
}
|
||||
|
||||
totali = nmax * m_numYlms;
|
||||
totali = (bigint) nmax * m_numYlms;
|
||||
for (bigint tn = 0; tn < totali; tn++) {
|
||||
m_clist_r[tn] = m_clist_r[tn] * double(weight);
|
||||
m_clist_i[tn] = m_clist_i[tn] * double(weight);
|
||||
@ -1174,9 +1173,10 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
|
||||
for (int n = 1; n < nmax + 1; n++) {
|
||||
int i = 0;
|
||||
for (int l = 0; l < lmax + 1; l++) {
|
||||
r_int = m_rip_array[(idpair - 1) * m_nmax * (m_lmax + 1) + (n - 1) * (m_lmax + 1) + l];
|
||||
r_int_temp =
|
||||
m_rip_darray[(idpair - 1) * m_nmax * (m_lmax + 1) + (n - 1) * (m_lmax + 1) + l];
|
||||
r_int = m_rip_array[(idpair - 1) * m_nmax * (m_lmax + 1) +
|
||||
(bigint) (n - 1) * (m_lmax + 1) + l];
|
||||
r_int_temp = m_rip_darray[(idpair - 1) * m_nmax * (m_lmax + 1) +
|
||||
(bigint) (n - 1) * (m_lmax + 1) + l];
|
||||
|
||||
for (int ii = 0; ii < 3; ii++) dr_int[ii] = r_int_temp * 2.0 * alpha * rvec[ii] / r;
|
||||
|
||||
|
||||
@ -81,8 +81,8 @@ PairRANN::PairRANN(LAMMPS *lmp) : Pair(lmp)
|
||||
manybody_flag = 1;
|
||||
allocated = 0;
|
||||
nelements = -1;
|
||||
elements = NULL;
|
||||
mass = NULL;
|
||||
elements = nullptr;
|
||||
mass = nullptr;
|
||||
|
||||
// set comm size needed by this Pair
|
||||
// comm unused for now.
|
||||
@ -304,7 +304,7 @@ void PairRANN::read_file(char *filename)
|
||||
while (eof == 0) {
|
||||
ptr=fgets(linetemp,longline,fp);
|
||||
linenum++;
|
||||
if (ptr == NULL) {
|
||||
if (ptr == nullptr) {
|
||||
if (check_potential()) {
|
||||
error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
|
||||
}
|
||||
@ -345,7 +345,7 @@ void PairRANN::read_file(char *filename)
|
||||
strtemp=utils::trim_comment(linetemp);
|
||||
linenum++;
|
||||
}
|
||||
if (ptr == NULL) {
|
||||
if (ptr == nullptr) {
|
||||
if (check_potential()) {
|
||||
error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
|
||||
}
|
||||
@ -532,7 +532,7 @@ void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string
|
||||
(*linenum)++;
|
||||
Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
|
||||
line1 = values1.as_vector();
|
||||
if (ptr==NULL)error->one(filename,*linenum,"unexpected end of potential file!");
|
||||
if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
|
||||
nwords = line1.size();
|
||||
if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
|
||||
for (k=0;k<net[l].dimensions[i];k++) {
|
||||
@ -575,7 +575,7 @@ void PairRANN::read_activation_functions(std::vector<std::string> line,std::vect
|
||||
for (l=0;l<nelements;l++) {
|
||||
if (line[1].compare(elements[l])==0) {
|
||||
if (net[l].layers==0)error->one(filename,linenum-1,"networklayers must be defined before activation functions.");
|
||||
i = strtol(line[2].c_str(),NULL,10);
|
||||
i = strtol(line[2].c_str(),nullptr,10);
|
||||
if (i>=net[l].layers || i<0)error->one(filename,linenum-1,"invalid activation layer");
|
||||
delete activation[l][i];
|
||||
activation[l][i]=create_activation(line1[0].c_str());
|
||||
@ -1232,10 +1232,8 @@ RANN::Fingerprint *PairRANN::create_fingerprint(const char *style)
|
||||
else if (strcmp(style,"bondspin")==0) {
|
||||
return new RANN::Fingerprint_bondspin(this);
|
||||
}
|
||||
char str[128];
|
||||
sprintf(str,"Unknown fingerprint style %s",style);
|
||||
error->one(FLERR,str);
|
||||
return NULL;
|
||||
error->one(FLERR,"Unknown fingerprint style {}",style);
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
@ -1247,9 +1245,7 @@ RANN::Activation *PairRANN::create_activation(const char *style)
|
||||
else if (strcmp(style,"sigI")==0) {
|
||||
return new RANN::Activation_sigI(this);
|
||||
}
|
||||
char str[128];
|
||||
sprintf(str,"Unknown activation style %s",style);
|
||||
error->one(FLERR,str);
|
||||
return NULL;
|
||||
error->one(FLERR,"Unknown activation style {}",style);
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
@ -70,6 +70,7 @@ PairSNAP::~PairSNAP()
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(scale);
|
||||
}
|
||||
|
||||
}
|
||||
@ -164,6 +165,7 @@ void PairSNAP::compute(int eflag, int vflag)
|
||||
// for neighbors of I within cutoff:
|
||||
// compute Fij = dEi/dRj = -dEi/dRi
|
||||
// add to Fi, subtract from Fj
|
||||
// scaling is that for type I
|
||||
|
||||
snaptr->compute_yi(beta[ii]);
|
||||
|
||||
@ -178,12 +180,12 @@ void PairSNAP::compute(int eflag, int vflag)
|
||||
|
||||
snaptr->compute_deidrj(fij);
|
||||
|
||||
f[i][0] += fij[0];
|
||||
f[i][1] += fij[1];
|
||||
f[i][2] += fij[2];
|
||||
f[j][0] -= fij[0];
|
||||
f[j][1] -= fij[1];
|
||||
f[j][2] -= fij[2];
|
||||
f[i][0] += fij[0]*scale[itype][itype];
|
||||
f[i][1] += fij[1]*scale[itype][itype];
|
||||
f[i][2] += fij[2]*scale[itype][itype];
|
||||
f[j][0] -= fij[0]*scale[itype][itype];
|
||||
f[j][1] -= fij[1]*scale[itype][itype];
|
||||
f[j][2] -= fij[2]*scale[itype][itype];
|
||||
|
||||
// tally per-atom virial contribution
|
||||
|
||||
@ -224,6 +226,7 @@ void PairSNAP::compute(int eflag, int vflag)
|
||||
}
|
||||
}
|
||||
}
|
||||
evdwl *= scale[itype][itype];
|
||||
ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
|
||||
}
|
||||
|
||||
@ -351,9 +354,9 @@ void PairSNAP::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag,n+1,n+1,"pair:setflag");
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
memory->create(scale,n+1,n+1,"pair:scale");
|
||||
map = new int[n+1];
|
||||
}
|
||||
|
||||
@ -408,11 +411,16 @@ void PairSNAP::coeff(int narg, char **arg)
|
||||
}
|
||||
|
||||
// Calculate maximum cutoff for all elements
|
||||
|
||||
rcutmax = 0.0;
|
||||
for (int ielem = 0; ielem < nelements; ielem++)
|
||||
rcutmax = MAX(2.0*radelem[ielem]*rcutfac,rcutmax);
|
||||
|
||||
// set default scaling
|
||||
int n = atom->ntypes;
|
||||
for (int ii = 0; ii < n+1; ii++)
|
||||
for (int jj = 0; jj < n+1; jj++)
|
||||
scale[ii][jj] = 1.0;
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -441,6 +449,7 @@ void PairSNAP::init_style()
|
||||
double PairSNAP::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
scale[j][i] = scale[i][j];
|
||||
return (radelem[map[i]] +
|
||||
radelem[map[j]])*rcutfac;
|
||||
}
|
||||
@ -711,6 +720,7 @@ double PairSNAP::memory_usage()
|
||||
int n = atom->ntypes+1;
|
||||
bytes += (double)n*n*sizeof(int); // setflag
|
||||
bytes += (double)n*n*sizeof(double); // cutsq
|
||||
bytes += (double)n*n*sizeof(double); // scale
|
||||
bytes += (double)n*sizeof(int); // map
|
||||
bytes += (double)beta_max*ncoeff*sizeof(double); // bispectrum
|
||||
bytes += (double)beta_max*ncoeff*sizeof(double); // beta
|
||||
@ -720,3 +730,9 @@ double PairSNAP::memory_usage()
|
||||
return bytes;
|
||||
}
|
||||
|
||||
void *PairSNAP::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"scale") == 0) return (void *) scale;
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
@ -34,6 +34,7 @@ class PairSNAP : public Pair {
|
||||
virtual void init_style();
|
||||
virtual double init_one(int, int);
|
||||
virtual double memory_usage();
|
||||
virtual void *extract(const char *, int &);
|
||||
|
||||
double rcutfac, quadraticflag; // declared public to workaround gcc 4.9
|
||||
int ncoeff; // compiler bug, manifest in KOKKOS package
|
||||
@ -55,6 +56,7 @@ class PairSNAP : public Pair {
|
||||
double **coeffelem; // element bispectrum coefficients
|
||||
double **beta; // betas for all atoms in list
|
||||
double **bispectrum; // bispectrum components for all atoms in list
|
||||
double **scale; // for thermodynamic integration
|
||||
int twojmax, switchflag, bzeroflag, bnormflag;
|
||||
int chemflag, wselfallflag;
|
||||
int chunksize;
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -16,24 +15,25 @@
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "omp_compat.h"
|
||||
#include "angle_fourier_simple_omp.h"
|
||||
#include <cmath>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include "omp_compat.h"
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define SMALL 0.001
|
||||
#define SMALL 0.0001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleFourierSimpleOMP::AngleFourierSimpleOMP(class LAMMPS *lmp)
|
||||
: AngleFourierSimple(lmp), ThrOMP(lmp,THR_ANGLE)
|
||||
AngleFourierSimpleOMP::AngleFourierSimpleOMP(class LAMMPS *lmp) :
|
||||
AngleFourierSimple(lmp), ThrOMP(lmp, THR_ANGLE)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
}
|
||||
@ -42,14 +42,14 @@ AngleFourierSimpleOMP::AngleFourierSimpleOMP(class LAMMPS *lmp)
|
||||
|
||||
void AngleFourierSimpleOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
const int nthreads = comm->nthreads;
|
||||
const int inum = neighbor->nanglelist;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
@ -62,34 +62,40 @@ void AngleFourierSimpleOMP::compute(int eflag, int vflag)
|
||||
if (inum > 0) {
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr);
|
||||
else eval<1,1,0>(ifrom, ito, thr);
|
||||
if (force->newton_bond)
|
||||
eval<1, 1, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 1, 0>(ifrom, ito, thr);
|
||||
} else {
|
||||
if (force->newton_bond) eval<1,0,1>(ifrom, ito, thr);
|
||||
else eval<1,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_bond)
|
||||
eval<1, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
} else {
|
||||
if (force->newton_bond) eval<0,0,1>(ifrom, ito, thr);
|
||||
else eval<0,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_bond)
|
||||
eval<0, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<0, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
}
|
||||
thr->timer(Timer::BOND);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_BOND>
|
||||
void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr)
|
||||
void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData *const thr)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double term,sgn;
|
||||
double rsq1,rsq2,r1,r2,c,cn,th,nth,a,a11,a12,a22;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double term, sgn;
|
||||
double rsq1, rsq2, r1, r2, c, cn, th, nth, a, a11, a12, a22;
|
||||
|
||||
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0];
|
||||
const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const int4_t *_noalias const anglelist = (int4_t *) neighbor->anglelist[0];
|
||||
const int nlocal = atom->nlocal;
|
||||
eangle = 0.0;
|
||||
|
||||
@ -105,7 +111,7 @@ void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr)
|
||||
dely1 = x[i1].y - x[i2].y;
|
||||
delz1 = x[i1].z - x[i2].z;
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -114,13 +120,13 @@ void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr)
|
||||
dely2 = x[i3].y - x[i2].y;
|
||||
delz2 = x[i3].z - x[i2].z;
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
@ -128,38 +134,38 @@ void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr)
|
||||
// force & energy
|
||||
|
||||
th = acos(c);
|
||||
nth = N[type]*acos(c);
|
||||
nth = N[type] * acos(c);
|
||||
cn = cos(nth);
|
||||
term = k[type]*(1.0+C[type]*cn);
|
||||
term = k[type] * (1.0 + C[type] * cn);
|
||||
|
||||
if (EFLAG) eangle = term;
|
||||
|
||||
// handle sin(n th)/sin(th) singulatiries
|
||||
|
||||
if (fabs(c)-1.0 > 0.0001) {
|
||||
a = k[type]*C[type]*N[type]*sin(nth)/sin(th);
|
||||
if (fabs(c) - 1.0 > SMALL) {
|
||||
a = k[type] * C[type] * N[type] * sin(nth) / sin(th);
|
||||
} else {
|
||||
if (c >= 0.0) {
|
||||
term = 1.0 - c;
|
||||
sgn = 1.0;
|
||||
} else {
|
||||
term = 1.0 + c;
|
||||
sgn = ( fmodf((float)(N[type]),2.0) == 0.0f )?-1.0:1.0;
|
||||
sgn = (fmod(N[type], 2.0) == 0.0) ? -1.0 : 1.0;
|
||||
}
|
||||
a = N[type]+N[type]*(1.0-N[type]*N[type])*term/3.0;
|
||||
a = k[type]*C[type]*N[type]*(double)(sgn)*a;
|
||||
a = N[type] + N[type] * (1.0 - N[type] * N[type]) * term / 3.0;
|
||||
a = k[type] * C[type] * N[type] * sgn * a;
|
||||
}
|
||||
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -181,7 +187,8 @@ void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr)
|
||||
f[i3].z += f3[2];
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this,i1,i2,i3,nlocal,NEWTON_BOND,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2,thr);
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this, i1, i2, i3, nlocal, NEWTON_BOND, eangle, f1, f3, delx1, dely1, delz1,
|
||||
delx2, dely2, delz2, thr);
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -19,21 +18,21 @@
|
||||
|
||||
#include "angle_fourier_simple.h"
|
||||
|
||||
#include <cmath>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define SMALL 0.001
|
||||
#define SMALL 0.0001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -60,14 +59,14 @@ AngleFourierSimple::~AngleFourierSimple()
|
||||
|
||||
void AngleFourierSimple::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double term,sgn;
|
||||
double rsq1,rsq2,r1,r2,c,cn,th,nth,a,a11,a12,a22;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double term, sgn;
|
||||
double rsq1, rsq2, r1, r2, c, cn, th, nth, a, a11, a12, a22;
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -88,7 +87,7 @@ void AngleFourierSimple::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -97,13 +96,13 @@ void AngleFourierSimple::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
@ -111,38 +110,38 @@ void AngleFourierSimple::compute(int eflag, int vflag)
|
||||
// force & energy
|
||||
|
||||
th = acos(c);
|
||||
nth = N[type]*acos(c);
|
||||
nth = N[type] * acos(c);
|
||||
cn = cos(nth);
|
||||
term = k[type]*(1.0+C[type]*cn);
|
||||
term = k[type] * (1.0 + C[type] * cn);
|
||||
|
||||
if (eflag) eangle = term;
|
||||
|
||||
// handle sin(n th)/sin(th) singulatiries
|
||||
|
||||
if (fabs(c)-1.0 > 0.0001) {
|
||||
a = k[type]*C[type]*N[type]*sin(nth)/sin(th);
|
||||
if (fabs(c) - 1.0 > SMALL) {
|
||||
a = k[type] * C[type] * N[type] * sin(nth) / sin(th);
|
||||
} else {
|
||||
if (c >= 0.0) {
|
||||
term = 1.0 - c;
|
||||
sgn = 1.0;
|
||||
} else {
|
||||
term = 1.0 + c;
|
||||
sgn = ( fmodf((float)(N[type]),2.0) == 0.0f )?-1.0:1.0;
|
||||
sgn = (fmod(N[type], 2.0) == 0.0) ? -1.0 : 1.0;
|
||||
}
|
||||
a = N[type]+N[type]*(1.0-N[type]*N[type])*term/3.0;
|
||||
a = k[type]*C[type]*N[type]*(double)(sgn)*a;
|
||||
a = N[type] + N[type] * (1.0 - N[type] * N[type]) * term / 3.0;
|
||||
a = k[type] * C[type] * N[type] * sgn * a;
|
||||
}
|
||||
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -164,8 +163,9 @@ void AngleFourierSimple::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -176,11 +176,11 @@ void AngleFourierSimple::allocate()
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(C,n+1,"angle:C");
|
||||
memory->create(N,n+1,"angle:N");
|
||||
memory->create(k, n + 1, "angle:k");
|
||||
memory->create(C, n + 1, "angle:C");
|
||||
memory->create(N, n + 1, "angle:N");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
memory->create(setflag, n + 1, "angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
@ -190,15 +190,15 @@ void AngleFourierSimple::allocate()
|
||||
|
||||
void AngleFourierSimple::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 4) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (narg != 4) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
double k_one = utils::numeric(FLERR,arg[1],false,lmp);
|
||||
double C_one = utils::numeric(FLERR,arg[2],false,lmp);
|
||||
double N_one = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
double C_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
double N_one = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
@ -209,14 +209,14 @@ void AngleFourierSimple::coeff(int narg, char **arg)
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleFourierSimple::equilibrium_angle(int i)
|
||||
{
|
||||
return (MY_PI/N[i]);
|
||||
return (MY_PI / N[i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -225,9 +225,9 @@ double AngleFourierSimple::equilibrium_angle(int i)
|
||||
|
||||
void AngleFourierSimple::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&C[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&N[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&C[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&N[1], sizeof(double), atom->nangletypes, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -239,13 +239,13 @@ void AngleFourierSimple::read_restart(FILE *fp)
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&C[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&N[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &C[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &N[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&C[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&N[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&C[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&N[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
@ -256,8 +256,7 @@ void AngleFourierSimple::read_restart(FILE *fp)
|
||||
|
||||
void AngleFourierSimple::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->nangletypes; i++)
|
||||
fprintf(fp,"%d %g %g %g\n",i,k[i],C[i],N[i]);
|
||||
for (int i = 1; i <= atom->nangletypes; i++) fprintf(fp, "%d %g %g %g\n", i, k[i], C[i], N[i]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -269,21 +268,21 @@ double AngleFourierSimple::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
double cn = cos(N[type]*acos(c));
|
||||
double cn = cos(N[type] * acos(c));
|
||||
|
||||
double eng = k[type]*(1.0+C[type]*cn);
|
||||
double eng = k[type] * (1.0 + C[type] * cn);
|
||||
return eng;
|
||||
}
|
||||
|
||||
@ -15,6 +15,7 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing Author: Xipeng Wang, Simon Ramirez-Hinestrosa
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_wf_cut.h"
|
||||
|
||||
#include "atom.h"
|
||||
@ -391,5 +392,5 @@ void *PairWFCut::extract(const char *str, int &dim)
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"nu") == 0) return (void *) nu;
|
||||
if (strcmp(str,"mu") == 0) return (void *) mu;
|
||||
return NULL;
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -13,32 +12,34 @@
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "compute_property_local.h"
|
||||
#include <cstring>
|
||||
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "update.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
|
||||
enum{TYPE,RADIUS};
|
||||
enum { NONE, NEIGH, PAIR, BOND, ANGLE, DIHEDRAL, IMPROPER };
|
||||
enum { TYPE, RADIUS };
|
||||
|
||||
#define DELTA 10000
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg),
|
||||
vlocal(NULL), alocal(NULL), indices(NULL), pack_choice(NULL)
|
||||
Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr), indices(nullptr),
|
||||
pack_choice(nullptr)
|
||||
{
|
||||
if (narg < 4) error->all(FLERR,"Illegal compute property/local command");
|
||||
if (narg < 4) error->all(FLERR, "Illegal compute property/local command");
|
||||
|
||||
local_flag = 1;
|
||||
nvalues = narg - 3;
|
||||
@ -50,220 +51,198 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
|
||||
nvalues = 0;
|
||||
int iarg = 3;
|
||||
while (iarg < narg) {
|
||||
i = iarg-3;
|
||||
i = iarg - 3;
|
||||
|
||||
if (strcmp(arg[iarg],"natom1") == 0) {
|
||||
if (strcmp(arg[iarg], "natom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_patom1;
|
||||
if (kindflag != NONE && kindflag != NEIGH)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = NEIGH;
|
||||
} else if (strcmp(arg[iarg],"natom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "natom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_patom2;
|
||||
if (kindflag != NONE && kindflag != NEIGH)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = NEIGH;
|
||||
} else if (strcmp(arg[iarg],"ntype1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "ntype1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_ptype1;
|
||||
if (kindflag != NONE && kindflag != NEIGH)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = NEIGH;
|
||||
} else if (strcmp(arg[iarg],"ntype2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "ntype2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_ptype2;
|
||||
if (kindflag != NONE && kindflag != NEIGH)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = NEIGH;
|
||||
|
||||
} else if (strcmp(arg[iarg],"patom1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "patom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_patom1;
|
||||
if (kindflag != NONE && kindflag != PAIR)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = PAIR;
|
||||
} else if (strcmp(arg[iarg],"patom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "patom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_patom2;
|
||||
if (kindflag != NONE && kindflag != PAIR)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = PAIR;
|
||||
} else if (strcmp(arg[iarg],"ptype1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "ptype1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_ptype1;
|
||||
if (kindflag != NONE && kindflag != PAIR)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = PAIR;
|
||||
} else if (strcmp(arg[iarg],"ptype2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "ptype2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_ptype2;
|
||||
if (kindflag != NONE && kindflag != PAIR)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = PAIR;
|
||||
|
||||
} else if (strcmp(arg[iarg],"batom1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "batom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_batom1;
|
||||
if (kindflag != NONE && kindflag != BOND)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = BOND;
|
||||
} else if (strcmp(arg[iarg],"batom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "batom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_batom2;
|
||||
if (kindflag != NONE && kindflag != BOND)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = BOND;
|
||||
} else if (strcmp(arg[iarg],"btype") == 0) {
|
||||
} else if (strcmp(arg[iarg], "btype") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_btype;
|
||||
if (kindflag != NONE && kindflag != BOND)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = BOND;
|
||||
|
||||
} else if (strcmp(arg[iarg],"aatom1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "aatom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_aatom1;
|
||||
if (kindflag != NONE && kindflag != ANGLE)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = ANGLE;
|
||||
} else if (strcmp(arg[iarg],"aatom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "aatom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_aatom2;
|
||||
if (kindflag != NONE && kindflag != ANGLE)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = ANGLE;
|
||||
} else if (strcmp(arg[iarg],"aatom3") == 0) {
|
||||
} else if (strcmp(arg[iarg], "aatom3") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_aatom3;
|
||||
if (kindflag != NONE && kindflag != ANGLE)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = ANGLE;
|
||||
} else if (strcmp(arg[iarg],"atype") == 0) {
|
||||
} else if (strcmp(arg[iarg], "atype") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_atype;
|
||||
if (kindflag != NONE && kindflag != ANGLE)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = ANGLE;
|
||||
|
||||
} else if (strcmp(arg[iarg],"datom1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "datom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_datom1;
|
||||
if (kindflag != NONE && kindflag != DIHEDRAL)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = DIHEDRAL;
|
||||
} else if (strcmp(arg[iarg],"datom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "datom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_datom2;
|
||||
if (kindflag != NONE && kindflag != DIHEDRAL)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = DIHEDRAL;
|
||||
} else if (strcmp(arg[iarg],"datom3") == 0) {
|
||||
} else if (strcmp(arg[iarg], "datom3") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_datom3;
|
||||
if (kindflag != NONE && kindflag != DIHEDRAL)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = DIHEDRAL;
|
||||
} else if (strcmp(arg[iarg],"datom4") == 0) {
|
||||
} else if (strcmp(arg[iarg], "datom4") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_datom4;
|
||||
if (kindflag != NONE && kindflag != DIHEDRAL)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = DIHEDRAL;
|
||||
} else if (strcmp(arg[iarg],"dtype") == 0) {
|
||||
} else if (strcmp(arg[iarg], "dtype") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_dtype;
|
||||
if (kindflag != NONE && kindflag != DIHEDRAL)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = DIHEDRAL;
|
||||
|
||||
} else if (strcmp(arg[iarg],"iatom1") == 0) {
|
||||
} else if (strcmp(arg[iarg], "iatom1") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_iatom1;
|
||||
if (kindflag != NONE && kindflag != IMPROPER)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = IMPROPER;
|
||||
} else if (strcmp(arg[iarg],"iatom2") == 0) {
|
||||
} else if (strcmp(arg[iarg], "iatom2") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_iatom2;
|
||||
if (kindflag != NONE && kindflag != IMPROPER)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = IMPROPER;
|
||||
} else if (strcmp(arg[iarg],"iatom3") == 0) {
|
||||
} else if (strcmp(arg[iarg], "iatom3") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_iatom3;
|
||||
if (kindflag != NONE && kindflag != IMPROPER)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = IMPROPER;
|
||||
} else if (strcmp(arg[iarg],"iatom4") == 0) {
|
||||
} else if (strcmp(arg[iarg], "iatom4") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_iatom4;
|
||||
if (kindflag != NONE && kindflag != IMPROPER)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = IMPROPER;
|
||||
} else if (strcmp(arg[iarg],"itype") == 0) {
|
||||
} else if (strcmp(arg[iarg], "itype") == 0) {
|
||||
pack_choice[i] = &ComputePropertyLocal::pack_itype;
|
||||
if (kindflag != NONE && kindflag != IMPROPER)
|
||||
error->all(FLERR,
|
||||
"Compute property/local cannot use these inputs together");
|
||||
error->all(FLERR, "Compute property/local cannot use these inputs together");
|
||||
kindflag = IMPROPER;
|
||||
|
||||
} else break;
|
||||
} else
|
||||
break;
|
||||
|
||||
nvalues++;
|
||||
iarg++;
|
||||
}
|
||||
|
||||
if (nvalues == 1) size_local_cols = 0;
|
||||
else size_local_cols = nvalues;
|
||||
if (nvalues == 1)
|
||||
size_local_cols = 0;
|
||||
else
|
||||
size_local_cols = nvalues;
|
||||
|
||||
// optional args
|
||||
|
||||
cutstyle = TYPE;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"cutoff") == 0) {
|
||||
if (iarg+2 > narg)
|
||||
error->all(FLERR,"Illegal compute property/local command");
|
||||
if (strcmp(arg[iarg+1],"type") == 0) cutstyle = TYPE;
|
||||
else if (strcmp(arg[iarg+1],"radius") == 0) cutstyle = RADIUS;
|
||||
else error->all(FLERR,"Illegal compute property/local command");
|
||||
if (strcmp(arg[iarg], "cutoff") == 0) {
|
||||
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute property/local command");
|
||||
if (strcmp(arg[iarg + 1], "type") == 0)
|
||||
cutstyle = TYPE;
|
||||
else if (strcmp(arg[iarg + 1], "radius") == 0)
|
||||
cutstyle = RADIUS;
|
||||
else
|
||||
error->all(FLERR, "Illegal compute property/local command");
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal compute property/local command");
|
||||
} else
|
||||
error->all(FLERR, "Illegal compute property/local command");
|
||||
}
|
||||
|
||||
// error check
|
||||
|
||||
if (atom->molecular == 2 && (kindflag == BOND || kindflag == ANGLE ||
|
||||
kindflag == DIHEDRAL || kindflag == IMPROPER))
|
||||
error->all(FLERR,"Compute property/local does not (yet) work "
|
||||
if (atom->molecular == 2 &&
|
||||
(kindflag == BOND || kindflag == ANGLE || kindflag == DIHEDRAL || kindflag == IMPROPER))
|
||||
error->all(FLERR,
|
||||
"Compute property/local does not (yet) work "
|
||||
"with atom_style template");
|
||||
|
||||
if (kindflag == BOND && atom->avec->bonds_allow == 0)
|
||||
error->all(FLERR,
|
||||
"Compute property/local for property that isn't allocated");
|
||||
error->all(FLERR, "Compute property/local for property that isn't allocated");
|
||||
if (kindflag == ANGLE && atom->avec->angles_allow == 0)
|
||||
error->all(FLERR,
|
||||
"Compute property/local for property that isn't allocated");
|
||||
error->all(FLERR, "Compute property/local for property that isn't allocated");
|
||||
if (kindflag == DIHEDRAL && atom->avec->dihedrals_allow == 0)
|
||||
error->all(FLERR,
|
||||
"Compute property/local for property that isn't allocated");
|
||||
error->all(FLERR, "Compute property/local for property that isn't allocated");
|
||||
if (kindflag == IMPROPER && atom->avec->impropers_allow == 0)
|
||||
error->all(FLERR,
|
||||
"Compute property/local for property that isn't allocated");
|
||||
error->all(FLERR, "Compute property/local for property that isn't allocated");
|
||||
if (cutstyle == RADIUS && !atom->radius_flag)
|
||||
error->all(FLERR,"Compute property/local requires atom attribute radius");
|
||||
error->all(FLERR, "Compute property/local requires atom attribute radius");
|
||||
|
||||
nmax = 0;
|
||||
vlocal = NULL;
|
||||
alocal = NULL;
|
||||
vlocal = nullptr;
|
||||
alocal = nullptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputePropertyLocal::~ComputePropertyLocal()
|
||||
{
|
||||
delete [] pack_choice;
|
||||
delete[] pack_choice;
|
||||
memory->destroy(vlocal);
|
||||
memory->destroy(alocal);
|
||||
memory->destroy(indices);
|
||||
@ -274,10 +253,10 @@ ComputePropertyLocal::~ComputePropertyLocal()
|
||||
void ComputePropertyLocal::init()
|
||||
{
|
||||
if (kindflag == NEIGH || kindflag == PAIR) {
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"No pair style is defined for compute property/local");
|
||||
if (force->pair == nullptr)
|
||||
error->all(FLERR, "No pair style is defined for compute property/local");
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute property/local");
|
||||
error->all(FLERR, "Pair style does not support compute property/local");
|
||||
}
|
||||
|
||||
// for NEIGH/PAIR need an occasional half neighbor list
|
||||
@ -285,7 +264,7 @@ void ComputePropertyLocal::init()
|
||||
// this should enable it to always be a copy list (e.g. for granular pstyle)
|
||||
|
||||
if (kindflag == NEIGH || kindflag == PAIR) {
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
@ -296,12 +275,18 @@ void ComputePropertyLocal::init()
|
||||
// do initial memory allocation so that memory_usage() is correct
|
||||
// cannot be done yet for NEIGH/PAIR, since neigh list does not exist
|
||||
|
||||
if (kindflag == NEIGH) ncount = 0;
|
||||
else if (kindflag == PAIR) ncount = 0;
|
||||
else if (kindflag == BOND) ncount = count_bonds(0);
|
||||
else if (kindflag == ANGLE) ncount = count_angles(0);
|
||||
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
|
||||
else if (kindflag == IMPROPER) ncount = count_impropers(0);
|
||||
if (kindflag == NEIGH)
|
||||
ncount = 0;
|
||||
else if (kindflag == PAIR)
|
||||
ncount = 0;
|
||||
else if (kindflag == BOND)
|
||||
ncount = count_bonds(0);
|
||||
else if (kindflag == ANGLE)
|
||||
ncount = count_angles(0);
|
||||
else if (kindflag == DIHEDRAL)
|
||||
ncount = count_dihedrals(0);
|
||||
else if (kindflag == IMPROPER)
|
||||
ncount = count_impropers(0);
|
||||
|
||||
if (ncount > nmax) reallocate(ncount);
|
||||
size_local_rows = ncount;
|
||||
@ -322,22 +307,34 @@ void ComputePropertyLocal::compute_local()
|
||||
|
||||
// count local entries and generate list of indices
|
||||
|
||||
if (kindflag == NEIGH) ncount = count_pairs(0,0);
|
||||
else if (kindflag == PAIR) ncount = count_pairs(0,1);
|
||||
else if (kindflag == BOND) ncount = count_bonds(0);
|
||||
else if (kindflag == ANGLE) ncount = count_angles(0);
|
||||
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
|
||||
else if (kindflag == IMPROPER) ncount = count_impropers(0);
|
||||
if (kindflag == NEIGH)
|
||||
ncount = count_pairs(0, 0);
|
||||
else if (kindflag == PAIR)
|
||||
ncount = count_pairs(0, 1);
|
||||
else if (kindflag == BOND)
|
||||
ncount = count_bonds(0);
|
||||
else if (kindflag == ANGLE)
|
||||
ncount = count_angles(0);
|
||||
else if (kindflag == DIHEDRAL)
|
||||
ncount = count_dihedrals(0);
|
||||
else if (kindflag == IMPROPER)
|
||||
ncount = count_impropers(0);
|
||||
|
||||
if (ncount > nmax) reallocate(ncount);
|
||||
size_local_rows = ncount;
|
||||
|
||||
if (kindflag == NEIGH) ncount = count_pairs(1,0);
|
||||
else if (kindflag == PAIR) ncount = count_pairs(1,1);
|
||||
else if (kindflag == BOND) ncount = count_bonds(1);
|
||||
else if (kindflag == ANGLE) ncount = count_angles(1);
|
||||
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(1);
|
||||
else if (kindflag == IMPROPER) ncount = count_impropers(1);
|
||||
if (kindflag == NEIGH)
|
||||
ncount = count_pairs(1, 0);
|
||||
else if (kindflag == PAIR)
|
||||
ncount = count_pairs(1, 1);
|
||||
else if (kindflag == BOND)
|
||||
ncount = count_bonds(1);
|
||||
else if (kindflag == ANGLE)
|
||||
ncount = count_angles(1);
|
||||
else if (kindflag == DIHEDRAL)
|
||||
ncount = count_dihedrals(1);
|
||||
else if (kindflag == IMPROPER)
|
||||
ncount = count_impropers(1);
|
||||
|
||||
// fill vector or array with local values
|
||||
|
||||
@ -346,8 +343,7 @@ void ComputePropertyLocal::compute_local()
|
||||
(this->*pack_choice[0])(0);
|
||||
} else {
|
||||
if (alocal) buf = &alocal[0][0];
|
||||
for (int n = 0; n < nvalues; n++)
|
||||
(this->*pack_choice[n])(n);
|
||||
for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n);
|
||||
}
|
||||
}
|
||||
|
||||
@ -361,10 +357,10 @@ void ComputePropertyLocal::compute_local()
|
||||
|
||||
int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
|
||||
{
|
||||
int i,j,m,ii,jj,inum,jnum,itype,jtype;
|
||||
tagint itag,jtag;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,radsum;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, m, ii, jj, inum, jnum, itype, jtype;
|
||||
tagint itag, jtag;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, rsq, radsum;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
@ -415,9 +411,9 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
|
||||
if (newton_pair == 0 && j >= nlocal) {
|
||||
jtag = tag[j];
|
||||
if (itag > jtag) {
|
||||
if ((itag+jtag) % 2 == 0) continue;
|
||||
if ((itag + jtag) % 2 == 0) continue;
|
||||
} else if (itag < jtag) {
|
||||
if ((itag+jtag) % 2 == 1) continue;
|
||||
if ((itag + jtag) % 2 == 1) continue;
|
||||
} else {
|
||||
if (x[j][2] < ztmp) continue;
|
||||
if (x[j][2] == ztmp) {
|
||||
@ -430,14 +426,14 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
if (forceflag) {
|
||||
if (cutstyle == TYPE) {
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
} else {
|
||||
radsum = radius[i] + radius[j];
|
||||
if (rsq >= radsum*radsum) continue;
|
||||
if (rsq >= radsum * radsum) continue;
|
||||
}
|
||||
}
|
||||
|
||||
@ -463,7 +459,7 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
|
||||
|
||||
int ComputePropertyLocal::count_bonds(int flag)
|
||||
{
|
||||
int i,atom1,atom2;
|
||||
int i, atom1, atom2;
|
||||
|
||||
int *num_bond = atom->num_bond;
|
||||
tagint **bond_atom = atom->bond_atom;
|
||||
@ -504,7 +500,7 @@ int ComputePropertyLocal::count_bonds(int flag)
|
||||
|
||||
int ComputePropertyLocal::count_angles(int flag)
|
||||
{
|
||||
int i,atom1,atom2,atom3;
|
||||
int i, atom1, atom2, atom3;
|
||||
|
||||
int *num_angle = atom->num_angle;
|
||||
tagint **angle_atom1 = atom->angle_atom1;
|
||||
@ -546,7 +542,7 @@ int ComputePropertyLocal::count_angles(int flag)
|
||||
|
||||
int ComputePropertyLocal::count_dihedrals(int flag)
|
||||
{
|
||||
int i,atom1,atom2,atom3,atom4;
|
||||
int i, atom1, atom2, atom3, atom4;
|
||||
|
||||
int *num_dihedral = atom->num_dihedral;
|
||||
tagint **dihedral_atom1 = atom->dihedral_atom1;
|
||||
@ -589,7 +585,7 @@ int ComputePropertyLocal::count_dihedrals(int flag)
|
||||
|
||||
int ComputePropertyLocal::count_impropers(int flag)
|
||||
{
|
||||
int i,atom1,atom2,atom3,atom4;
|
||||
int i, atom1, atom2, atom3, atom4;
|
||||
|
||||
int *num_improper = atom->num_improper;
|
||||
tagint **improper_atom1 = atom->improper_atom1;
|
||||
@ -633,16 +629,16 @@ void ComputePropertyLocal::reallocate(int n)
|
||||
|
||||
if (nvalues == 1) {
|
||||
memory->destroy(vlocal);
|
||||
memory->create(vlocal,nmax,"property/local:vector_local");
|
||||
memory->create(vlocal, nmax, "property/local:vector_local");
|
||||
vector_local = vlocal;
|
||||
} else {
|
||||
memory->destroy(alocal);
|
||||
memory->create(alocal,nmax,nvalues,"property/local:array_local");
|
||||
memory->create(alocal, nmax, nvalues, "property/local:array_local");
|
||||
array_local = alocal;
|
||||
}
|
||||
|
||||
memory->destroy(indices);
|
||||
memory->create(indices,nmax,2,"property/local:indices");
|
||||
memory->create(indices, nmax, 2, "property/local:indices");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -651,8 +647,8 @@ void ComputePropertyLocal::reallocate(int n)
|
||||
|
||||
double ComputePropertyLocal::memory_usage()
|
||||
{
|
||||
double bytes = (double)nmax*nvalues * sizeof(double);
|
||||
bytes += (double)nmax*2 * sizeof(int);
|
||||
double bytes = (double) nmax * nvalues * sizeof(double);
|
||||
bytes += (double) nmax * 2 * sizeof(int);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
@ -736,7 +732,7 @@ void ComputePropertyLocal::pack_batom1(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_batom2(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **bond_atom = atom->bond_atom;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -751,7 +747,7 @@ void ComputePropertyLocal::pack_batom2(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_btype(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
int **bond_type = atom->bond_type;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -766,7 +762,7 @@ void ComputePropertyLocal::pack_btype(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_aatom1(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **angle_atom1 = atom->angle_atom1;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -781,7 +777,7 @@ void ComputePropertyLocal::pack_aatom1(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_aatom2(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **angle_atom2 = atom->angle_atom2;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -796,7 +792,7 @@ void ComputePropertyLocal::pack_aatom2(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_aatom3(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **angle_atom3 = atom->angle_atom3;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -811,7 +807,7 @@ void ComputePropertyLocal::pack_aatom3(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_atype(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
int **angle_type = atom->angle_type;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -826,7 +822,7 @@ void ComputePropertyLocal::pack_atype(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_datom1(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **dihedral_atom1 = atom->dihedral_atom1;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -841,7 +837,7 @@ void ComputePropertyLocal::pack_datom1(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_datom2(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **dihedral_atom2 = atom->dihedral_atom2;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -856,7 +852,7 @@ void ComputePropertyLocal::pack_datom2(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_datom3(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **dihedral_atom3 = atom->dihedral_atom3;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -871,7 +867,7 @@ void ComputePropertyLocal::pack_datom3(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_datom4(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **dihedral_atom4 = atom->dihedral_atom4;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -886,7 +882,7 @@ void ComputePropertyLocal::pack_datom4(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_dtype(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
int **dihedral_type = atom->dihedral_type;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -901,7 +897,7 @@ void ComputePropertyLocal::pack_dtype(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_iatom1(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **improper_atom1 = atom->improper_atom1;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -916,7 +912,7 @@ void ComputePropertyLocal::pack_iatom1(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_iatom2(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **improper_atom2 = atom->improper_atom2;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -931,7 +927,7 @@ void ComputePropertyLocal::pack_iatom2(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_iatom3(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **improper_atom3 = atom->improper_atom3;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -946,7 +942,7 @@ void ComputePropertyLocal::pack_iatom3(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_iatom4(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
tagint **improper_atom4 = atom->improper_atom4;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
@ -961,7 +957,7 @@ void ComputePropertyLocal::pack_iatom4(int n)
|
||||
|
||||
void ComputePropertyLocal::pack_itype(int n)
|
||||
{
|
||||
int i,j;
|
||||
int i, j;
|
||||
int **improper_type = atom->improper_type;
|
||||
|
||||
for (int m = 0; m < ncount; m++) {
|
||||
|
||||
Reference in New Issue
Block a user