tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2886 from akohlmey/pair-coul-exclude

Browse Source 
Add a pair style coul/exclude for removing excluded coulomb interactions
This commit is contained in:
Axel Kohlmeyer
2021-08-19 17:25:57 -04:00
committed by GitHub
parent ee5273fbb1 8aeaa82909
commit 537f0eed8e
38 changed files with 1050 additions and 561 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_pair.rst
View File

@ -75,6 +75,7 @@ OPT.
* :doc:`coul/debye (gko) <pair_coul>`
* :doc:`coul/diel (o) <pair_coul_diel>`
* :doc:`coul/dsf (gko) <pair_coul>`
* :doc:`coul/exclude <pair_coul>`
* :doc:`coul/long (gko) <pair_coul>`
* :doc:`coul/long/cs (g) <pair_cs>`
* :doc:`coul/long/dielectric <pair_dielectric>`

2
doc/src/Developer_utils.rst
View File

@ -336,7 +336,7 @@ repetitive tasks.
The :cpp:class:`LAMMPS_NS::ArgInfo` class provides an abstraction
for parsing references to compute or fix styles, variables or custom
integer or double properties handled by :doc:`fix property/atom <fix_property_atom>`.
These would start with a "c\_", "f\_", "v\_", "d\_", "d2\_", "i\_", or "i2_"
These would start with a "c\_", "f\_", "v\_", "d\_", "d2\_", "i\_", or "i2\_"
followed by the ID or name of than instance and may be postfixed with
one or two array indices "[<number>]" with numbers > 0.

21
doc/src/pair_coul.rst
View File

@ -10,6 +10,7 @@
.. index:: pair_style coul/dsf/gpu
.. index:: pair_style coul/dsf/kk
.. index:: pair_style coul/dsf/omp
.. index:: pair_style coul/exclude
.. index:: pair_style coul/cut/global
.. index:: pair_style coul/cut/global/omp
.. index:: pair_style coul/long
@ -42,6 +43,9 @@ pair_style coul/dsf command
Accelerator Variants: *coul/dsf/gpu*, *coul/dsf/kk*, *coul/dsf/omp*
pair_style coul/exclude command
===============================
pair_style coul/cut/global command
==================================
@ -83,6 +87,7 @@ Syntax
pair_style coul/cut cutoff
pair_style coul/debye kappa cutoff
pair_style coul/dsf alpha cutoff
pair_style coul/exclude cutoff
pair_style coul/cut/global cutoff
pair_style coul/long cutoff
pair_style coul/wolf alpha cutoff
@ -110,6 +115,9 @@ Examples
pair_style coul/dsf 0.05 10.0
pair_coeff * *
pair_style hybrid/overlay coul/exclude 10.0 ...
pair_coeff * * coul/exclude
pair_style coul/long 10.0
pair_coeff * *
@ -257,6 +265,19 @@ as style *coul/cut* except that it allows only a single global cutoff
and thus makes it compatible for use in combination with long-range
coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
Pair style *coul/exclude* computes Coulombic interactions like *coul/cut*
but **only** applies them to excluded pairs using a scaling factor
of :math:`\gamma - 1.0` with :math:`\gamma` being the factor assigned
to that excluded pair via the :doc:`special_bonds coul <special_bonds>`
setting. With this it is possible to treat Coulomb interactions for
molecular systems with :doc:`kspace style scafacos <kspace_style>`,
which always computes the *full* Coulomb interactions without exclusions.
Pair style *coul/exclude* will then *subtract* the excluded interactions
accordingly. So to achieve the same forces as with ``pair_style lj/cut/coul/long 12.0``
with ``kspace_style pppm 1.0e-6``, one would use
``pair_style hybrid/overlay lj/cut 12.0 coul/exclude 12.0`` with
``kspace_style scafacos p3m 1.0e-6``.
Styles *coul/long* and *coul/msm* compute the same Coulombic
interactions as style *coul/cut* except that an additional damping
factor is applied so it can be used in conjunction with the

1
doc/src/pair_style.rst
View File

@ -139,6 +139,7 @@ accelerated styles exist.
* :doc:`coul/debye <pair_coul>` - cutoff Coulomb potential with Debye screening
* :doc:`coul/diel <pair_coul_diel>` - Coulomb potential with dielectric permittivity
* :doc:`coul/dsf <pair_coul>` - Coulomb with damped-shifted-force model
* :doc:`coul/exclude <pair_coul>` - subtract Coulomb potential for excluded pairs
* :doc:`coul/long <pair_coul>` - long-range Coulomb potential
* :doc:`coul/long/cs <pair_cs>` - long-range Coulomb potential and core/shell
* :doc:`coul/long/dielectric <pair_dielectric>` -

307
src/EXTRA-PAIR/pair_coul_exclude.cpp Normal file
View File

@ -0,0 +1,307 @@
// 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.
------------------------------------------------------------------------- */
#include "pair_coul_exclude.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCoulExclude::PairCoulExclude(LAMMPS *lmp) : Pair(lmp) {
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairCoulExclude::~PairCoulExclude()
{
if (copymode) return;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
void PairCoulExclude::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double rsq,r2inv,rinv,forcecoul,factor_coul;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0;
ev_init(eflag,vflag);
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
// skip over all non-excluded pairs and subtract excluded coulomb
if (sbmask(j) == 0) continue;
factor_coul = special_coul[sbmask(j)] - 1.0;
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
forcecoul = qqrd2e * qtmp*q[j]*rinv;
fpair = factor_coul*forcecoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag)
ecoul = factor_coul * qqrd2e * qtmp*q[j]*rinv;
if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0,ecoul,fpair,delx,dely,delz);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairCoulExclude::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCoulExclude::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_global = utils::numeric(FLERR,arg[0],false,lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairCoulExclude::coeff(int narg, char **arg)
{
if (narg != 2)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCoulExclude::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style coul/exclude requires atom attribute q");
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCoulExclude::init_one(int i, int j)
{
return cut_global;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulExclude::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulExclude::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) {
utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulExclude::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulExclude::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairCoulExclude::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d\n",i);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairCoulExclude::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d\n",i,j);
}
/* ---------------------------------------------------------------------- */
double PairCoulExclude::single(int i, int j, int /*itype*/, int /*jtype*/,
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r2inv,rinv,forcecoul,phicoul;
factor_coul -= 1.0; // we want to subtract the excluded coulomb
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv;
fforce = factor_coul*forcecoul * r2inv;
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv;
return factor_coul*phicoul;
}
/* ---------------------------------------------------------------------- */
void *PairCoulExclude::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_global;
return nullptr;
}

72
src/EXTRA-PAIR/pair_coul_exclude.h Normal file
View File

@ -0,0 +1,72 @@
/* -*- 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 PAIR_CLASS
// clang-format off
PairStyle(coul/exclude,PairCoulExclude);
// clang-format on
#else
#ifndef LMP_PAIR_COUL_EXCLUDE_H
#define LMP_PAIR_COUL_EXCLUDE_H
#include "pair.h"
namespace LAMMPS_NS {
class PairCoulExclude : public Pair {
public:
PairCoulExclude(class LAMMPS *);
virtual ~PairCoulExclude();
virtual void compute(int, int);
virtual void settings(int, char **);
virtual void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual void write_data(FILE *);
virtual void write_data_all(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
virtual void *extract(const char *, int &);
protected:
double cut_global;
virtual void allocate();
};
} // 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.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style coul/cut requires atom attribute q
The atom style defined does not have these attributes.
*/

216
src/KSPACE/pair_coul_long.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -18,27 +17,27 @@
#include "pair_coul_long.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
@ -69,16 +68,16 @@ PairCoulLong::~PairCoulLong()
void PairCoulLong::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itable,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double fraction,table;
double r,r2inv,forcecoul,factor_coul;
double grij,expm2,prefactor,t,erfc;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itable, itype, jtype;
double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair;
double fraction, table;
double r, r2inv, forcecoul, factor_coul;
double grij, expm2, prefactor, t, erfc;
int *ilist, *jlist, *numneigh, **firstneigh;
double rsq;
ecoul = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -114,58 +113,57 @@ void PairCoulLong::compute(int eflag, int vflag)
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 (rsq < cut_coulsq) {
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
expm2 = exp(-grij * grij);
t = 1.0 / (1.0 + EWALD_P * grij);
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
prefactor = qqrd2e * scale[itype][jtype] * qtmp * q[j] / r;
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
table = ftable[itable] + fraction * dftable[itable];
forcecoul = scale[itype][jtype] * qtmp * q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
table = ctable[itable] + fraction * dctable[itable];
prefactor = scale[itype][jtype] * qtmp * q[j] * table;
forcecoul -= (1.0 - factor_coul) * prefactor;
}
}
fpair = forcecoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
if (!ncoultablebits || rsq <= tabinnersq)
ecoul = prefactor*erfc;
ecoul = prefactor * erfc;
else {
table = etable[itable] + fraction*detable[itable];
ecoul = scale[itype][jtype] * qtmp*q[j] * table;
table = etable[itable] + fraction * detable[itable];
ecoul = scale[itype][jtype] * qtmp * q[j] * table;
}
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz);
}
}
}
@ -182,14 +180,13 @@ void PairCoulLong::allocate()
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
for (int j = i; j <= n; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(scale,n+1,n+1,"pair:scale");
memory->create(scale, n + 1, n + 1, "pair:scale");
}
/* ----------------------------------------------------------------------
@ -198,9 +195,9 @@ void PairCoulLong::allocate()
void PairCoulLong::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
if (narg != 1) error->all(FLERR, "Illegal pair_style command");
cut_coul = utils::numeric(FLERR,arg[0],false,lmp);
cut_coul = utils::numeric(FLERR, arg[0], false, lmp);
}
/* ----------------------------------------------------------------------
@ -209,23 +206,23 @@ void PairCoulLong::settings(int narg, char **arg)
void PairCoulLong::coeff(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg != 2) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
scale[i][j] = 1.0;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -234,22 +231,20 @@ void PairCoulLong::coeff(int narg, char **arg)
void PairCoulLong::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q");
if (!atom->q_flag) error->all(FLERR, "Pair style lj/cut/coul/long requires atom attribute q");
neighbor->request(this,instance_me);
neighbor->request(this, instance_me);
cut_coulsq = cut_coul * cut_coul;
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == nullptr)
error->all(FLERR,"Pair style requires a KSpace style");
if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
// setup force tables
if (ncoultablebits) init_tables(cut_coul,nullptr);
if (ncoultablebits) init_tables(cut_coul, nullptr);
}
/* ----------------------------------------------------------------------
@ -259,7 +254,7 @@ void PairCoulLong::init_style()
double PairCoulLong::init_one(int i, int j)
{
scale[j][i] = scale[i][j];
return cut_coul+2.0*qdist;
return cut_coul + 2.0 * qdist;
}
/* ----------------------------------------------------------------------
@ -272,9 +267,8 @@ void PairCoulLong::write_restart(FILE *fp)
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j])
fwrite(&scale[i][j],sizeof(double),1,fp);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) fwrite(&scale[i][j], sizeof(double), 1, fp);
}
}
@ -288,15 +282,15 @@ void PairCoulLong::read_restart(FILE *fp)
allocate();
int i,j;
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) utils::sfread(FLERR,&scale[i][j],sizeof(double),1,fp,nullptr,error);
MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world);
if (me == 0) utils::sfread(FLERR, &scale[i][j], sizeof(double), 1, fp, nullptr, error);
MPI_Bcast(&scale[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
@ -307,11 +301,11 @@ void PairCoulLong::read_restart(FILE *fp)
void PairCoulLong::write_restart_settings(FILE *fp)
{
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&ncoultablebits,sizeof(int),1,fp);
fwrite(&tabinner,sizeof(double),1,fp);
fwrite(&cut_coul, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&ncoultablebits, sizeof(int), 1, fp);
fwrite(&tabinner, sizeof(double), 1, fp);
}
/* ----------------------------------------------------------------------
@ -321,63 +315,61 @@ void PairCoulLong::write_restart_settings(FILE *fp)
void PairCoulLong::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world);
MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world);
}
/* ---------------------------------------------------------------------- */
double PairCoulLong::single(int i, int j, int /*itype*/, int /*jtype*/,
double rsq,
double factor_coul, double /*factor_lj*/,
double &fforce)
double PairCoulLong::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq,
double factor_coul, double /*factor_lj*/, double &fforce)
{
double r2inv,r,grij,expm2,t,erfc,prefactor;
double fraction,table,forcecoul,phicoul;
double r2inv, r, grij, expm2, t, erfc, prefactor;
double fraction, table, forcecoul, phicoul;
int itable;
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
expm2 = exp(-grij * grij);
t = 1.0 / (1.0 + EWALD_P * grij);
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
table = ftable[itable] + fraction * dftable[itable];
forcecoul = atom->q[i] * atom->q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = atom->q[i]*atom->q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
table = ctable[itable] + fraction * dctable[itable];
prefactor = atom->q[i] * atom->q[j] * table;
forcecoul -= (1.0 - factor_coul) * prefactor;
}
}
fforce = forcecoul * r2inv;
if (!ncoultablebits || rsq <= tabinnersq)
phicoul = prefactor*erfc;
phicoul = prefactor * erfc;
else {
table = etable[itable] + fraction*detable[itable];
phicoul = atom->q[i]*atom->q[j] * table;
table = etable[itable] + fraction * detable[itable];
phicoul = atom->q[i] * atom->q[j] * table;
}
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
return phicoul;
}
@ -386,11 +378,11 @@ double PairCoulLong::single(int i, int j, int /*itype*/, int /*jtype*/,
void *PairCoulLong::extract(const char *str, int &dim)
{
if (strcmp(str,"cut_coul") == 0) {
if (strcmp(str, "cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"scale") == 0) {
if (strcmp(str, "scale") == 0) {
dim = 2;
return (void *) scale;
}

6
src/ML-RANN/rann_fingerprint_bond.cpp
View File

@ -591,8 +591,8 @@ void Fingerprint_bond::do3bodyfeatureset_doubleneighborloop(double * features,do
delz = zn[jj];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq>rc*rc) {
expr[jj][0]=0;
continue;
expr[jj][0]=0;
continue;
}
double r1 = (rsq*((double)res)*cutinv2);
int m1 = (int)r1;
@ -603,7 +603,7 @@ void Fingerprint_bond::do3bodyfeatureset_doubleneighborloop(double * features,do
double *p2 = &expcuttable[(m1+1)*kmax];
double *p3 = &expcuttable[(m1+2)*kmax];
for (kk=0;kk<kmax;kk++) {
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
}
double* q = &dfctable[m1-1];
double* r2 = &rinvsqrttable[m1-1];

12
src/ML-RANN/rann_fingerprint_bondscreened.cpp
View File

@ -627,16 +627,16 @@ void Fingerprint_bondscreened::do3bodyfeatureset_doubleneighborloop(double * fea
if (Bij[jj]==false) {continue;}
jtype = tn[jj];
if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
expr[jj][0]=0;
continue;
expr[jj][0]=0;
continue;
}
delx = xn[jj];
dely = yn[jj];
delz = zn[jj];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq>rc*rc) {
expr[jj][0]=0;
continue;
expr[jj][0]=0;
continue;
}
double r1 = (rsq*((double)res)*cutinv2);
int m1 = (int)r1;
@ -647,8 +647,8 @@ void Fingerprint_bondscreened::do3bodyfeatureset_doubleneighborloop(double * fea
double *p2 = &expcuttable[(m1+1)*kmax];
double *p3 = &expcuttable[(m1+2)*kmax];
for (kk=0;kk<kmax;kk++) {
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
expr[jj][kk] *= Sik[jj];
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
expr[jj][kk] *= Sik[jj];
}
double* q = &dfctable[m1-1];
double* r2 = &rinvsqrttable[m1-1];

50
src/OPENMP/npair_full_multi_omp.cpp
View File

@ -100,7 +100,7 @@ void NPairFullMultiOmp::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
@ -110,34 +110,34 @@ void NPairFullMultiOmp::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (i == j) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (i == j) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}

16
src/OPENMP/npair_half_multi_newtoff_omp.cpp
View File

@ -102,7 +102,7 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -114,16 +114,16 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {

110
src/OPENMP/npair_half_multi_newton_omp.cpp
View File

@ -101,7 +101,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// if same size: uses half stencil so check central bin
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -117,41 +117,41 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
for (j = js; j >= 0; j = bins[j]) {
for (j = js; j >= 0; j = bins[j]) {
if ((icollection != jcollection) && (j < i)) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
// for all collections, loop over all atoms in other bins in stencil, store every pair
@ -159,38 +159,38 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
// stencil is half if i same size as j
// stencil is full if i smaller than j
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[i] = i;

56
src/OPENMP/npair_half_multi_newton_tri_omp.cpp
View File

@ -102,7 +102,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -113,12 +113,12 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -133,30 +133,30 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[i] = i;

36
src/OPENMP/npair_half_size_multi_newtoff_omp.cpp
View File

@ -92,7 +92,7 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -104,27 +104,27 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >=0; j = bins[j]) {
if (j <= i) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >=0; j = bins[j]) {
if (j <= i) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}

80
src/OPENMP/npair_half_size_multi_newton_omp.cpp
View File

@ -91,7 +91,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// if same size: uses half stencil so check central bin
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -107,33 +107,33 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
for (j = js; j >= 0; j = bins[j]) {
for (j = js; j >= 0; j = bins[j]) {
if(icollection != jcollection and j < i) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
@ -142,31 +142,31 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
// stencil is half if i same size as j
// stencil is full if i smaller than j
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}
ilist[i] = i;

40
src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp
View File

@ -92,7 +92,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
@ -104,12 +104,12 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -126,21 +126,21 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}
ilist[i] = i;

23
src/SCAFACOS/scafacos.cpp
View File

@ -23,6 +23,7 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "pair_hybrid.h"
#include <cstdlib>
#include <cstring>
@ -108,10 +109,26 @@ void Scafacos::init()
if (domain->triclinic) error->all(FLERR, "Cannot use ScaFaCoS with triclinic domain yet");
if (atom->natoms > INT_MAX && sizeof(int) != 8)
error->all(FLERR, "Scafacos atom count exceeds 2B");
error->all(FLERR, "ScaFaCoS atom count exceeds 2B");
if (atom->molecular != Atom::ATOMIC)
error->all(FLERR, "Cannot use Scafacos with molecular charged systems yet");
if ((atom->molecular != Atom::ATOMIC) && (atom->nbonds + atom->nangles + atom->ndihedrals) > 0) {
int flag = 0;
if ((force->special_coul[1] == 1.0) && (force->special_coul[2] == 1.0) &&
(force->special_coul[3] == 1.0))
++flag;
PairHybrid *ph = reinterpret_cast<PairHybrid *>(force->pair_match("^hybrid", 0));
if (ph) {
for (int isub = 0; isub < ph->nstyles; ++isub)
if (force->pair_match("coul/exclude", 0, isub)) ++flag;
} else {
if (force->pair_match("coul/exclude", 0)) ++flag;
}
if (!flag)
error->all(FLERR,
"Must use pair style coul/exclude or 'special_bonds coul 1.0 1.0 1.0'"
"for molecular charged systems with ScaFaCos");
}
FCSResult result;

4
src/nbin_multi.cpp
View File

@ -202,7 +202,7 @@ void NBinMulti::setup_bins(int /*style*/)
// test for too many global bins in any dimension due to huge global domain
if (bbox[0]*binsizeinv > MAXSMALLINT || bbox[1]*binsizeinv > MAXSMALLINT ||
bbox[2]*binsizeinv > MAXSMALLINT)
bbox[2]*binsizeinv > MAXSMALLINT)
error->all(FLERR,"Domain too large for neighbor bins");
// create actual bins
@ -233,7 +233,7 @@ void NBinMulti::setup_bins(int /*style*/)
bininvz_multi[n] = 1.0 / binsizez_multi[n];
if (binsize_optimal*bininvx_multi[n] > CUT2BIN_RATIO ||
binsize_optimal*bininvy_multi[n] > CUT2BIN_RATIO)
binsize_optimal*bininvy_multi[n] > CUT2BIN_RATIO)
error->all(FLERR,"Cannot use neighbor bins - box size << cutoff");
if ((dimension == 3) && (binsize_optimal*bininvz_multi[n] > CUT2BIN_RATIO))
error->all(FLERR,"Cannot use neighbor bins - box size << cutoff");

2
src/neighbor.cpp
View File

@ -528,6 +528,7 @@ void Neighbor::init()
for (int isub=0; isub < ph->nstyles; ++isub) {
if (force->pair_match("coul/wolf",0,isub)
|| force->pair_match("coul/dsf",0,isub)
|| force->pair_match("coul/exclude",0)
|| force->pair_match("thole",0,isub))
++flag;
}
@ -536,6 +537,7 @@ void Neighbor::init()
} else {
if (force->pair_match("coul/wolf",0)
|| force->pair_match("coul/dsf",0)
|| force->pair_match("coul/exclude",0)
|| force->pair_match("thole",0))
special_flag[1] = special_flag[2] = special_flag[3] = 2;
}

48
src/npair_full_multi.cpp
View File

@ -89,7 +89,7 @@ void NPairFullMulti::build(NeighList *list)
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
@ -99,34 +99,34 @@ void NPairFullMulti::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (i == j) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (i == j) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}

14
src/npair_half_multi_newtoff.cpp
View File

@ -91,7 +91,7 @@ void NPairHalfMultiNewtoff::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -103,16 +103,16 @@ void NPairHalfMultiNewtoff::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {

110
src/npair_half_multi_newton.cpp
View File

@ -90,7 +90,7 @@ void NPairHalfMultiNewton::build(NeighList *list)
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// if same size: uses half stencil so check central bin
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -106,41 +106,41 @@ void NPairHalfMultiNewton::build(NeighList *list)
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
for (j = js; j >= 0; j = bins[j]) {
for (j = js; j >= 0; j = bins[j]) {
if((icollection != jcollection) && (j < i)) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
// for all collections, loop over all atoms in other bins in stencil, store every pair
@ -148,38 +148,38 @@ void NPairHalfMultiNewton::build(NeighList *list)
// stencil is half if i same size as j
// stencil is full if i smaller than j
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[inum++] = i;

54
src/npair_half_multi_newton_tri.cpp
View File

@ -90,7 +90,7 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -101,12 +101,12 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -123,28 +123,28 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[inum++] = i;

34
src/npair_half_size_multi_newtoff.cpp
View File

@ -81,7 +81,7 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -93,27 +93,27 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list)
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
if (j <= i) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}

76
src/npair_half_size_multi_newton.cpp
View File

@ -78,7 +78,7 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// if same size: uses half stencil so check central bin
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -94,33 +94,33 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
for (j = js; j >= 0; j = bins[j]) {
for (j = js; j >= 0; j = bins[j]) {
if ((icollection != jcollection) && (j < i)) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
@ -129,30 +129,30 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
// stencil is half if i same size as j
// stencil is full if i smaller than j
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}

40
src/npair_half_size_multi_newton_tri.cpp
View File

@ -78,7 +78,7 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -89,12 +89,12 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
@ -111,21 +111,21 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
if (rsq <= cutdistsq) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
neighptr[n++] = j;
}
}
}
}
ilist[inum++] = i;

4
src/nstencil.cpp
View File

@ -335,7 +335,7 @@ void NStencil::create_setup()
for (i = 0; i < n; ++i) {
stencil_multi[i] = new int*[n]();
for (j = 0; j < n; ++j) {
maxstencil_multi[i][j] = 0;
maxstencil_multi[i][j] = 0;
nstencil_multi[i][j] = 0;
stencil_multi[i][j] = nullptr;
}
@ -390,7 +390,7 @@ void NStencil::create_setup()
if(stencil_multi[i][j])
memory->destroy(stencil_multi[i][j]);
memory->create(stencil_multi[i][j], smax,
"neighstencil::stencil_multi");
"neighstencil::stencil_multi");
}
}
}

2
src/nstencil_full_multi_2d.cpp
View File

@ -77,7 +77,7 @@ void NStencilFullMulti2d::create()
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
nstencil_multi[icollection][jcollection] = ns;
}

4
src/nstencil_full_multi_3d.cpp
View File

@ -80,8 +80,8 @@ void NStencilFullMulti3d::create()
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
k*mbiny*mbinx + j*mbinx + i;
nstencil_multi[icollection][jcollection] = ns;

4
src/nstencil_half_multi_2d.cpp
View File

@ -92,12 +92,12 @@ void NStencilHalfMulti2d::create()
if (j > 0 || (j == 0 && i > 0)) {
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
}
}
} else {
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
}
nstencil_multi[icollection][jcollection] = ns;

6
src/nstencil_half_multi_2d_tri.cpp
View File

@ -90,12 +90,12 @@ void NStencilHalfMulti2dTri::create()
for (j = 0; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
} else {
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
if (bin_distance_multi(i,j,0,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j*mbinx + i;
}
nstencil_multi[icollection][jcollection] = ns;

12
src/nstencil_half_multi_3d.cpp
View File

@ -92,17 +92,17 @@ void NStencilHalfMulti3d::create()
for (k = 0; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (k > 0 || j > 0 || (j == 0 && i > 0)) {
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
if (k > 0 || j > 0 || (j == 0 && i > 0)) {
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
k*mbiny*mbinx + j*mbinx + i;
}
}
} else {
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
k*mbiny*mbinx + j*mbinx + i;
}

6
src/nstencil_half_multi_3d_tri.cpp
View File

@ -93,14 +93,14 @@ void NStencilHalfMulti3dTri::create()
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
stencil_multi[icollection][jcollection][ns++] =
k*mbiny*mbinx + j*mbinx + i;
} else {
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] =
k*mbiny*mbinx + j*mbinx + i;
}

4
src/pair.h
View File

@ -143,8 +143,6 @@ class Pair : protected Pointers {
void v_tally2_newton(int, double *, double *);
void v_tally3(int, int, int, double *, double *, double *, double *);
void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *);
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double,
double);
// general child-class methods
@ -274,6 +272,8 @@ class Pair : protected Pointers {
void ev_tally4(int, int, int, int, double, double *, double *, double *, double *, double *,
double *);
void ev_tally_tip4p(int, int *, double *, double, double);
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double,
double);
void v_tally2(int, int, double, double *);
void v_tally_tensor(int, int, int, int, double, double, double, double, double, double);
void virial_fdotr_compute();

11
src/pair_hybrid.cpp
View File

@ -581,18 +581,17 @@ void PairHybrid::init_style()
for (i = 1; i < 4; ++i) {
if (((force->special_lj[i] == 0.0) || (force->special_lj[i] == 1.0))
&& (force->special_lj[i] != special_lj[istyle][i]))
error->all(FLERR,"Pair_modify special setting for pair hybrid "
"incompatible with global special_bonds setting");
error->all(FLERR,"Pair_modify special lj 1-{} setting for pair hybrid substyle {} "
"incompatible with global special_bonds setting", i+1, keywords[istyle]);
}
}
if (special_coul[istyle]) {
for (i = 1; i < 4; ++i) {
if (((force->special_coul[i] == 0.0)
|| (force->special_coul[i] == 1.0))
if (((force->special_coul[i] == 0.0) || (force->special_coul[i] == 1.0))
&& (force->special_coul[i] != special_coul[istyle][i]))
error->all(FLERR,"Pair_modify special setting for pair hybrid "
"incompatible with global special_bonds setting");
error->all(FLERR,"Pair_modify special coul 1-{} setting for pair hybrid substyle {} "
"incompatible with global special_bonds setting", i+1, keywords[istyle]);
}
}
}

5
src/pair_hybrid.h
View File

@ -28,13 +28,14 @@ class PairHybrid : public Pair {
friend class ComputeSpin;
friend class FixGPU;
friend class FixIntel;
friend class FixOMP;
friend class FixNVESpin;
friend class FixOMP;
friend class Force;
friend class Info;
friend class Neighbor;
friend class Respa;
friend class PairDeprecated;
friend class Respa;
friend class Scafacos;
public:
PairHybrid(class LAMMPS *);

35
src/special.cpp
View File

@ -237,10 +237,8 @@ void Special::onetwo_build_newton()
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous), 0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous), (void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
@ -371,10 +369,8 @@ void Special::onethree_build()
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous), 0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous), (void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
@ -475,10 +471,8 @@ void Special::onefour_build()
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous), 0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous), (void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
@ -898,10 +892,8 @@ void Special::angle_trim()
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous), 0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous), (void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
@ -1111,10 +1103,8 @@ void Special::dihedral_trim()
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous), 0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous), (void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
@ -1246,9 +1236,8 @@ int Special::rendezvous_ids(int n, char *inbuf,
outbuf = same list of N PairRvous datums, routed to different procs
------------------------------------------------------------------------- */
int Special::rendezvous_pairs(int n, char *inbuf,
int &flag, int *&proclist, char *&outbuf,
void *ptr)
int Special::rendezvous_pairs(int n, char *inbuf, int &flag, int *&proclist,
char *&outbuf, void *ptr)
{
Special *sptr = (Special *) ptr;
Atom *atom = sptr->atom;

87
unittest/force-styles/tests/mol-pair-coul_exclude.yaml Normal file
View File

@ -0,0 +1,87 @@
---
lammps_version: 30 Jul 2021
date_generated: Thu Aug 19 05:44:56 2021
epsilon: 5e-14
skip_tests:
prerequisites: ! |
atom full
pair coul/exclude
pre_commands: ! ""
post_commands: ! |
pair_modify mix arithmetic
input_file: in.fourmol
pair_style: coul/exclude 8.0
pair_coeff: ! |
* *
extract: ! |
cut_coul 0
natoms: 29
init_vdwl: 0
init_coul: 976.7662140166423
init_stress: ! |2-
3.5388948320536065e+02 2.5335487380105781e+02 3.6952185701022364e+02 1.3686424356022002e+02 8.7895961622382046e+01 -5.5657699252762978e+00
init_forces: ! |2
1 -2.1370477255985712e+01 -1.4166598766893758e+01 2.9310440052146326e+01
2 2.3264732430469657e+01 1.5692001202883613e+01 -2.8641449819420849e+01
3 2.2449042305615680e-01 1.4668967160440745e+00 5.1210039450899436e-01
4 -7.0490395531893513e-01 -1.3264719224865684e+00 -8.1217625652904202e-01
5 -8.8255242947367551e-01 -2.5499496555354710e+00 -1.1751860565740935e-01
6 -2.5789282115617354e+01 2.5922083376147825e+01 3.9284668116535428e+01
7 6.5887213922755370e+00 -1.1379710951282684e+01 -4.1569698240745815e+01
8 1.0165793608296209e+01 -1.9732745647894593e+01 -3.1168434006878172e+01
9 6.2278611600582163e+00 7.5490085322004052e+00 3.6053658984732728e+01
10 2.9613956094579899e+00 -6.6704022762728634e+00 7.0753413211731742e-01
11 7.6925511893250798e-01 -1.0660633894890499e+00 1.5916232833331401e+00
12 -2.6730309790123474e+00 1.0151505797685552e+00 2.4428922002007205e+00
13 4.2028736668256519e+00 -2.3945364969381782e+00 2.2608110571641327e-01
14 -1.0444805762876637e+00 -2.1692594974964866e-01 -4.5361664712630878e+00
15 1.0294632170564109e+00 4.2802145197407322e+00 1.1966180837917026e+00
16 1.6139881756442762e+01 -1.2595226879230463e+01 -4.8159451500674258e+01
17 -1.9109741071175417e+01 1.6173277008988080e+01 4.3679278548085868e+01
18 -5.5233108291167312e+00 -3.1467038458483600e+01 1.1751414157603500e+02
19 -4.8039195406996384e+01 -2.8078667518868457e+01 -6.9999907578978252e+01
20 5.3562506236113109e+01 5.9545705977352057e+01 -4.7514233997056749e+01
21 -3.2000540242879481e+01 -3.9305792536444351e+01 1.1575486824768517e+02
22 -4.1655033980718414e+01 -6.9415565144852787e+00 -8.1937970834963011e+01
23 7.3655574223597895e+01 4.6247349050929628e+01 -3.3816897412722149e+01
24 2.5133287180289841e+01 -1.0428031336996401e+02 6.0768019640810614e+01
25 -6.6553531734748447e+01 2.0322561261478373e+01 -5.9750075573899906e+01
26 4.1420244554458598e+01 8.3957752108485636e+01 -1.0179440669107107e+00
27 1.6184779133783991e+01 -1.1478028188716922e+02 4.0330521475487721e+01
28 -6.9426784473130638e+01 4.0320458958275808e+01 -4.9831480196977843e+01
29 5.3242005339346647e+01 7.4459822928893402e+01 9.5009587214901217e+00
run_vdwl: 0
run_coul: 976.5089740760076
run_stress: ! |2-
3.5429311069548726e+02 2.5338890902230719e+02 3.6882695435821336e+02 1.3697994178745301e+02 8.8376897487416215e+01 -4.7604760003034263e+00
run_forces: ! |2
1 -2.1374527865789219e+01 -1.4202468442180830e+01 2.9210793601587412e+01
2 2.3263825454132370e+01 1.5734731460617787e+01 -2.8542033985914678e+01
3 2.2462585109904795e-01 1.4683529055153852e+00 5.1240421813750503e-01
4 -7.0182458619184773e-01 -1.3247801518720248e+00 -8.0930989810165044e-01
5 -8.8184812299538406e-01 -2.5499921086303359e+00 -1.1763132809557786e-01
6 -2.5769520757000304e+01 2.5925263738943752e+01 3.9219964964514880e+01
7 6.5906098601140588e+00 -1.1393221625638784e+01 -4.1479783467188192e+01
8 1.0136770975699770e+01 -1.9687788806065537e+01 -3.1157177158436330e+01
9 6.2325131117568873e+00 7.5029136682374347e+00 3.6011266563600110e+01
10 2.9627918492957752e+00 -6.6677419862633922e+00 7.1111494624978899e-01
11 7.6684242658501378e-01 -1.0650907890115915e+00 1.5865397504820036e+00
12 -2.6640875675365687e+00 1.0092746068515972e+00 2.4201156771591061e+00
13 4.2020911135059107e+00 -2.3836918937522000e+00 2.2652733116436718e-01
14 -1.0435191142174847e+00 -2.2501596985534378e-01 -4.5261002002201831e+00
15 1.0232065551846687e+00 4.2815515515559763e+00 1.2073921411042297e+00
16 1.6133578122629647e+01 -1.2611941628990227e+01 -4.8150714173103566e+01
17 -1.9101527306272342e+01 1.6189645470538316e+01 4.3676631017060792e+01
18 -5.0679097498234569e+00 -3.0925377041375377e+01 1.1693070865584599e+02
19 -4.8160088791916564e+01 -2.8252999531559283e+01 -6.9934716985727633e+01
20 5.3227998541740021e+01 5.9178376572934660e+01 -4.6995991670118372e+01
21 -3.2023145472732139e+01 -3.9029459539316420e+01 1.1553537552520324e+02
22 -4.1759300667456344e+01 -7.0799675491211342e+00 -8.1817190857460076e+01
23 7.3782446140188483e+01 4.6109427088437556e+01 -3.3718184667743166e+01
24 2.5304219247242145e+01 -1.0428415260557270e+02 6.0810156089628336e+01
25 -6.6772835122638725e+01 2.0222909012153103e+01 -5.9945279201780089e+01
26 4.1468615875396580e+01 8.4061243593419590e+01 -8.6487688784824801e-01
27 1.6328783769519674e+01 -1.1474592950014997e+02 4.0113097823544649e+01
28 -6.9456403419277834e+01 4.0310251785996094e+01 -4.9735673433675416e+01
29 5.3127619649758152e+01 7.4435677714153883e+01 9.6225756101307667e+00
...