ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2012-10-02 16:36:24 +00:00
parent 157d510873
commit 953061d658
18 changed files with 64 additions and 4389 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
src/KSPACE/Install.sh
View File

@ -3,17 +3,23 @@
if (test $1 = 1) then if (test $1 = 1) then
cp ewald.cpp .. cp ewald.cpp ..
cp ewald_n.cpp ..
cp pppm.cpp .. cp pppm.cpp ..
cp pppm_old.cpp .. cp pppm_old.cpp ..
cp pppm_cg.cpp .. cp pppm_cg.cpp ..
cp pppm_disp.cpp ..
cp pppm_disp_tip4p.cpp ..
cp pppm_tip4p.cpp .. cp pppm_tip4p.cpp ..
cp msm.cpp .. cp msm.cpp ..
cp pair_born_coul_long.cpp .. cp pair_born_coul_long.cpp ..
cp pair_buck_coul_long.cpp .. cp pair_buck_coul_long.cpp ..
cp pair_buck_coul.cpp ..
cp pair_coul_long.cpp .. cp pair_coul_long.cpp ..
cp pair_lj_cut_coul_long.cpp .. cp pair_lj_cut_coul_long.cpp ..
cp pair_lj_cut_coul_long_tip4p.cpp .. cp pair_lj_cut_coul_long_tip4p.cpp ..
cp pair_lj_charmm_coul_long.cpp .. cp pair_lj_charmm_coul_long.cpp ..
cp pair_lj_coul.cpp ..
cp pair_lj_coul_tip4p.cpp ..
cp pair_born_coul_msm.cpp .. cp pair_born_coul_msm.cpp ..
cp pair_buck_coul_msm.cpp .. cp pair_buck_coul_msm.cpp ..
cp pair_coul_msm.cpp .. cp pair_coul_msm.cpp ..
@ -25,18 +31,26 @@ if (test $1 = 1) then
cp remap_wrap.cpp .. cp remap_wrap.cpp ..
cp ewald.h .. cp ewald.h ..
cp ewald_n.h ..
cp kissfft.h .. cp kissfft.h ..
cp math_complex.h ..
cp math_vector.h ..
cp msm.h .. cp msm.h ..
cp pppm.h .. cp pppm.h ..
cp pppm_old.h .. cp pppm_old.h ..
cp pppm_cg.h .. cp pppm_cg.h ..
cp pppm_disp.h ..
cp pppm_disp_tip4p.h ..
cp pppm_tip4p.h .. cp pppm_tip4p.h ..
cp pair_born_coul_long.h .. cp pair_born_coul_long.h ..
cp pair_buck_coul_long.h .. cp pair_buck_coul_long.h ..
cp pair_buck_coul.h ..
cp pair_coul_long.h .. cp pair_coul_long.h ..
cp pair_lj_cut_coul_long.h .. cp pair_lj_cut_coul_long.h ..
cp pair_lj_cut_coul_long_tip4p.h .. cp pair_lj_cut_coul_long_tip4p.h ..
cp pair_lj_charmm_coul_long.h .. cp pair_lj_charmm_coul_long.h ..
cp pair_lj_coul.h ..
cp pair_lj_coul_tip4p.h ..
cp pair_born_coul_msm.h .. cp pair_born_coul_msm.h ..
cp pair_buck_coul_msm.h .. cp pair_buck_coul_msm.h ..
cp pair_coul_msm.h .. cp pair_coul_msm.h ..
@ -50,17 +64,23 @@ if (test $1 = 1) then
elif (test $1 = 0) then elif (test $1 = 0) then
rm -f ../ewald.cpp rm -f ../ewald.cpp
rm -f ../ewald_n.cpp
rm -f ../msm.cpp rm -f ../msm.cpp
rm -f ../pppm.cpp rm -f ../pppm.cpp
rm -f ../pppm_old.cpp rm -f ../pppm_old.cpp
rm -f ../pppm_cg.cpp rm -f ../pppm_cg.cpp
rm -f ../pppm_disp.cpp
rm -f ../pppm_disp_tip4p.cpp
rm -f ../pppm_tip4p.cpp rm -f ../pppm_tip4p.cpp
rm -f ../pair_born_coul_long.cpp rm -f ../pair_born_coul_long.cpp
rm -f ../pair_buck_coul_long.cpp rm -f ../pair_buck_coul_long.cpp
rm -f ../pair_buck_coul.cpp
rm -f ../pair_coul_long.cpp rm -f ../pair_coul_long.cpp
rm -f ../pair_lj_cut_coul_long.cpp rm -f ../pair_lj_cut_coul_long.cpp
rm -f ../pair_lj_cut_coul_long_tip4p.cpp rm -f ../pair_lj_cut_coul_long_tip4p.cpp
rm -f ../pair_lj_charmm_coul_long.cpp rm -f ../pair_lj_charmm_coul_long.cpp
rm -f ../pair_lj_coul.cpp
rm -f ../pair_lj_coul_tip4p.cpp
rm -f ../pair_born_coul_msm.cpp rm -f ../pair_born_coul_msm.cpp
rm -f ../pair_buck_coul_msm.cpp rm -f ../pair_buck_coul_msm.cpp
rm -f ../pair_coul_msm.cpp rm -f ../pair_coul_msm.cpp
@ -72,18 +92,26 @@ elif (test $1 = 0) then
rm -f ../remap_wrap.cpp rm -f ../remap_wrap.cpp
rm -f ../ewald.h rm -f ../ewald.h
rm -f ../ewald_n.h
rm -f ../kissfft.h rm -f ../kissfft.h
rm -f ../math_vector.h
rm -f ../math_complex.h
rm -f ../msm.h rm -f ../msm.h
rm -f ../pppm.h rm -f ../pppm.h
rm -f ../pppm_old.h rm -f ../pppm_old.h
rm -f ../pppm_cg.h rm -f ../pppm_cg.h
rm -f ../pppm_disp.h
rm -f ../pppm_disp_tip4p.h
rm -f ../pppm_tip4p.h rm -f ../pppm_tip4p.h
rm -f ../pair_born_coul_long.h rm -f ../pair_born_coul_long.h
rm -f ../pair_buck_coul_long.h rm -f ../pair_buck_coul_long.h
rm -f ../pair_buck_coul.h
rm -f ../pair_coul_long.h rm -f ../pair_coul_long.h
rm -f ../pair_lj_cut_coul_long.h rm -f ../pair_lj_cut_coul_long.h
rm -f ../pair_lj_cut_coul_long_tip4p.h rm -f ../pair_lj_cut_coul_long_tip4p.h
rm -f ../pair_lj_charmm_coul_long.h rm -f ../pair_lj_charmm_coul_long.h
rm -f ../pair_lj_coul.h
rm -f ../pair_lj_coul_tip4p.h
rm -f ../pair_born_coul_msm.h rm -f ../pair_born_coul_msm.h
rm -f ../pair_buck_coul_msm.h rm -f ../pair_buck_coul_msm.h
rm -f ../pair_coul_msm.h rm -f ../pair_coul_msm.h

3
src/Makefile
View File

@ -18,8 +18,7 @@ PACKAGE = asphere class2 colloid dipole fld gpu granular kim \
rigid shock srd xtc rigid shock srd xtc
PACKUSER = user-misc user-atc user-awpmd user-cg-cmm user-colvars \ PACKUSER = user-misc user-atc user-awpmd user-cg-cmm user-colvars \
user-cuda user-eff user-ewaldn user-omp user-molfile \ user-cuda user-eff user-omp user-molfile user-reaxc user-sph
user-reaxc user-sph
PACKLIB = gpu kim meam poems reax user-atc user-awpmd user-colvars \ PACKLIB = gpu kim meam poems reax user-atc user-awpmd user-colvars \
user-cuda user-molfile user-cuda user-molfile

33
src/USER-EWALDN/Install.sh
View File

@ -1,33 +0,0 @@
# Install/unInstall package files in LAMMPS
if (test $1 = 1) then
cp -p ewald_n.cpp ..
cp -p pair_buck_coul.cpp ..
#cp -p pair_lj_dipole.cpp ..
cp -p pair_lj_coul.cpp ..
cp -p ewald_n.h ..
cp -p pair_buck_coul.h ..
#cp -p pair_lj_dipole.h ..
cp -p pair_lj_coul.h ..
cp -p math_vector.h ..
cp -p math_complex.h ..
elif (test $1 = 0) then
rm -f ../ewald_n.cpp
rm -f ../pair_buck_coul.cpp
#rm -f ../pair_lj_dipole.cpp
rm -f ../pair_lj_coul.cpp
rm -f ../ewald_n.h
rm -f ../pair_buck_coul.h
#rm -f ../pair_lj_dipole.h
rm -f ../pair_lj_coul.h
rm -f ../math_vector.h
rm -f ../math_complex.h
fi

29
src/USER-EWALDN/README
View File

@ -1,29 +0,0 @@
This package implements 3 commands which can be used in a LAMMPS input
script: pair_style lj/coul, pair_style buck/coul, and kspace_style
ewald/n.
The "kspace_style ewald/n" command is similar to standard Ewald for
charges, but also enables the Lennard-Jones interaction, or any 1/r^N
interaction to be of infinite extent, instead of being cutoff. LAMMPS
pair potentials for long-range Coulombic interactions, such as
lj/cut/coul/long can be used with ewald/n. The two new pair_style
commands provide the modifications for the short-range LJ and
Buckingham interactions that can also be used with ewald/n.
Two other advantages of kspace_style ewald/n are that
a) it can be used with non-orthogonal (triclinic symmetry) simulation
boxes
b) it can include long-range summations not just for Coulombic
interactions (1/r), but also for dispersion interactions (1/r^6) and
dipole interactions (1/r^3).
Neither of these options is currently possible for other kspace styles
such as PPPM and ewald.
See the doc pages for these commands for details.
The person who created these files is Pieter in' t Veld while at
Sandia. He is now at BASF (pieter.intveld at basf.com). Contact him
directly if you have questions.

1221
src/USER-EWALDN/ewald_n.cpp
View File

File diff suppressed because it is too large Load Diff

95
src/USER-EWALDN/ewald_n.h
View File

@ -1,95 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef KSPACE_CLASS
KSpaceStyle(ewald/n,EwaldN)
#else
#ifndef LMP_EWALD_N_H
#define LMP_EWALD_N_H
#include "kspace.h"
#include "math_complex.h"
namespace LAMMPS_NS {
#define EWALD_NORDER 6
#define EWALD_NFUNCS 4
#define EWALD_MAX_NSUMS 10
#define EWALD_NSUMS {1, 1, 7, 1}
typedef struct cvector { complex x, y, z; } cvector;
typedef struct hvector { double x, y, z; } hvector;
typedef struct kvector { long x, y, z; } kvector;
class EwaldN : public KSpace {
public:
EwaldN(class LAMMPS *, int, char **);
~EwaldN();
void init();
void setup();
void compute(int, int);
double memory_usage() {return bytes;}
private:
double unit[6];
int function[EWALD_NFUNCS], first_output;
int nkvec, nkvec_max, nevec, nevec_max,
nbox, nfunctions, nsums, sums;
int peratom_allocate_flag;
int nmax;
double bytes;
double gsqmx,q2,b2;
double *kenergy, energy_self[EWALD_NFUNCS];
double *kvirial, virial_self[EWALD_NFUNCS];
double **energy_self_peratom;
double **virial_self_peratom;
cvector *ekr_local;
hvector *hvec;
kvector *kvec;
double mumurd2e, dielectric, *B, volume;
struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS];
complex *cek_local, *cek_global;
double rms(int, double, bigint, double, double);
void reallocate();
void allocate_peratom();
void reallocate_atoms();
void deallocate();
void deallocate_peratom();
void coefficients();
void init_coeffs();
void init_coeff_sums();
void init_self();
void init_self_peratom();
void compute_ek();
void compute_force();
void compute_surface();
void compute_energy();
void compute_energy_peratom();
void compute_virial();
void compute_virial_peratom();
void compute_slabcorr();
double NewtonSolve(double, double, bigint, double, double);
double f(double, double, bigint, double, double);
double derivf(double, double, bigint, double, double);
};
}
#endif
#endif

72
src/USER-EWALDN/math_complex.h
View File

@ -1,72 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
#ifndef LMP_MATH_COMPLEX_H
#define LMP_MATH_COMPLEX_H
#define COMPLEX_NULL {0, 0}
namespace LAMMPS_NS {
typedef struct complex {
double re, im; } complex;
}
#define C_MULT(d, x, y) { \
d.re = x.re*y.re-x.im*y.im; \
d.im = x.re*y.im+x.im*y.re; }
#define C_RMULT(d, x, y) { \
register complex t = x; \
d.re = t.re*y.re-t.im*y.im; \
d.im = t.re*y.im+t.im*y.re; }
#define C_CRMULT(d, x, y) { \
register complex t = x; \
d.re = t.re*y.re-t.im*y.im; \
d.im = -t.re*y.im-t.im*y.re; }
#define C_SMULT(d, x, y) { \
d.re = x.re*y; \
d.im = x.im*y; }
#define C_ADD(d, x, y) { \
d.re = x.re+y.re; \
d.im = x.im+y.im; }
#define C_SUBTR(d, x, y) { \
d.re = x.re-y.re; \
d.im = x.im-y.im; }
#define C_CONJ(d, x) { \
d.re = x.re; \
d.im = -x.im; }
#define C_SET(d, x, y) { \
d.re = x; \
d.im = y; }
#define C_ANGLE(d, angle) { \
register double a = angle; \
d.re = cos(a); \
d.im = sin(a); }
#define C_COPY(d, x) { \
memcpy(&d, &x, sizeof(complex)); }
#endif

448
src/USER-EWALDN/math_vector.h
View File

@ -1,448 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
#ifndef LMP_MATH_VECTOR_H
#define LMP_MATH_VECTOR_H
#include "math.h"
#include "string.h"
#define VECTOR_NULL {0, 0, 0}
#define SHAPE_NULL {0, 0, 0, 0, 0, 0}
#define FORM_NULL {0, 0, 0, 0, 0, 0}
#define MATRIX_NULL {VECTOR_NULL, VECTOR_NULL, VECTOR_NULL}
#define VECTOR4_NULL {0, 0, 0, 0}
#define QUATERNION_NULL {0, 0, 0, 0}
#define FORM4_NULL {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
#define FZERO 1e-15
#define fzero(x) (((x)>-FZERO) && ((x)<FZERO))
namespace LAMMPS_NS {
typedef double vector[3]; // 0:x 1:y 2:z
typedef int lvector[3];
typedef double shape[6]; // 0:xx 1:yy 2:zz 3:zy 4:zx 5:yx
typedef int lshape[6]; // xy=0 xz=0 yz=0;
typedef double form[6]; // 0:xx 1:yy 2:zz 3:zy 4:zx 5:yx
typedef int lform[6]; // xy=yx xz=zx yz=zy;
typedef vector matrix[3]; // 00:xx 11:yy 22:zz 21:zy 20:zx 10:yx
typedef lvector lmatrix[3];
typedef double vector4[4]; // 4D vector
typedef double form4[10]; // 0:00 1:11 2:22 3:33 4:32
// 5:31 6:30 7:21 8:20 9:10
// 01=10 02=20 03=30 etc
typedef double quaternion[4]; // quaternion
// vector operators
inline double vec_dot(vector &a, vector &b) { // a.b
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
inline void vec_null(vector &dest) {
memset(dest, 0, sizeof(vector)); }
inline void vec_neg(vector &dest) { // -a
dest[0] = -dest[0];
dest[1] = -dest[1];
dest[2] = -dest[2]; }
inline void vec_norm(vector &dest) { // a/|a|
register double f = sqrt(vec_dot(dest, dest));
dest[0] /= f;
dest[1] /= f;
dest[2] /= f; }
inline void vec_add(vector &dest, vector &src) { // a+b
dest[0] += src[0];
dest[1] += src[1];
dest[2] += src[2]; }
inline void vec_subtr(vector &dest, vector &src) { // a-b
dest[0] -= src[0];
dest[1] -= src[1];
dest[2] -= src[2]; }
inline void vec_mult(vector &dest, vector &src) { // a*b
dest[0] *= src[0];
dest[1] *= src[1];
dest[2] *= src[2]; }
inline void vec_div(vector &dest, vector &src) { // a/b
dest[0] /= src[0];
dest[1] /= src[1];
dest[2] /= src[2]; }
inline void vec_cross(vector &dest, vector &src) { // a x b
vector t = {
dest[1]*src[2]-dest[2]*src[1],
dest[2]*src[0]-dest[0]*src[2],
dest[0]*src[1]-dest[1]*src[0]};
memcpy(dest, t, sizeof(vector)); }
inline void vec_scalar_mult(vector &dest, double f) { // f*a
dest[0] *= f;
dest[1] *= f;
dest[2] *= f; }
inline void vec_to_lvec(lvector &dest, vector &src) {
dest[0] = (int) src[0];
dest[1] = (int) src[1];
dest[2] = (int) src[2]; }
inline void lvec_to_vec(vector &dest, lvector &src) {
dest[0] = (double) src[0];
dest[1] = (double) src[1];
dest[2] = (double) src[2]; }
// shape operators
inline double shape_det(shape &s) {
return s[0]*s[1]*s[2]; }
inline void shape_null(shape &dest) {
memset(dest, 0, sizeof(shape)); }
inline void shape_unit(shape &dest) {
memset(dest, 0, sizeof(shape));
dest[0] = dest[1] = dest[2] = 1.0; }
inline void shape_add(shape &dest, shape &src) { // h_a+h_b
dest[0] += src[0]; dest[1] += src[1]; dest[2] += src[2];
dest[3] += src[3]; dest[4] += src[4]; dest[5] += src[5]; }
inline void shape_subtr(shape &dest, shape &src) { // h_a-h_b
dest[0] -= src[0]; dest[1] -= src[1]; dest[2] -= src[2];
dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
inline void shape_inv(shape &h_inv, shape &h) { // h^-1
h_inv[0] = 1.0/h[0]; h_inv[1] = 1.0/h[1]; h_inv[2] = 1.0/h[2];
h_inv[3] = -h[3]/(h[1]*h[2]);
h_inv[4] = (h[3]*h[5]-h[1]*h[4])/(h[0]*h[1]*h[2]);
h_inv[5] = -h[5]/(h[0]*h[1]); }
inline void shape_dot(shape &dest, shape &src) { // h_a.h_b
dest[3] = dest[1]*src[3]+dest[3]*src[2];
dest[4] = dest[0]*src[4]+dest[5]*src[3]+dest[4]*src[2];
dest[5] = dest[0]*src[5]+dest[5]*src[1];
dest[0] *= src[0]; dest[1] *= src[1]; dest[2] *= src[2]; }
inline void shape_scalar_mult(shape &dest, double f) { // f*h
dest[0] *= f; dest[1] *= f; dest[2] *= f;
dest[3] *= f; dest[4] *= f; dest[5] *= f; }
inline void shape_vec_dot(vector &dest, vector &src, shape &h) {// h.a
dest[0] = h[0]*src[0]+h[5]*src[1]+h[4]*src[2];
dest[1] = h[1]*src[1]+h[3]*src[2];
dest[2] = h[2]*src[2]; }
inline void vec_shape_dot(vector &dest, shape &h, vector &src) {// a.h
dest[2] = h[4]*src[0]+h[3]*src[1]+h[2]*src[2];
dest[1] = h[5]*src[0]+h[1]*src[1];
dest[0] = h[0]*src[0]; }
inline void shape_to_matrix(matrix &dest, shape &h) { // m = h
dest[0][0] = h[0]; dest[1][0] = h[5]; dest[2][0] = h[4];
dest[0][1] = 0.0; dest[1][1] = h[1]; dest[2][1] = h[3];
dest[0][2] = 0.0; dest[1][2] = 0.0; dest[2][2] = h[2]; }
inline void shape_to_lshape(lshape &dest, shape &src) {
dest[0] = (long)src[0]; dest[1] = (long)src[1]; dest[2] = (long)src[2];
dest[3] = (long)src[3]; dest[4] = (long)src[4]; dest[5] = (long)src[5]; }
inline void lshape_to_shape(shape &dest, lshape &src) {
dest[0] = (double)src[0]; dest[1] = (double)src[1]; dest[2] = (double)src[2];
dest[3] = (double)src[3]; dest[4] = (double)src[4]; dest[5] = (double)src[5];}
// form operators
inline double form_det(form &m) { // |m|
return m[0]*(m[1]*m[2]-m[3]*m[3])+
m[4]*(2.0*m[3]*m[5]-m[1]*m[4])-m[2]*m[5]*m[5]; }
inline void form_null(form &dest) {
memset(dest, 0, sizeof(form)); }
inline void form_unit(form &dest) {
memset(dest, 0, sizeof(form));
dest[0] = dest[1] = dest[2] = 1.0; }
inline void form_add(form &dest, form &src) { // m_a+m_b
dest[0] += src[0]; dest[1] += src[1]; dest[2] += src[2];
dest[3] += src[3]; dest[4] += src[4]; dest[5] += src[5]; }
inline void form_subtr(shape &dest, form &src) { // m_a-m_b
dest[0] -= src[0]; dest[1] -= src[1]; dest[2] -= src[2];
dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
inline int form_inv(form &m_inv, form &m) { // m^-1
register double det = form_det(m);
if (fzero(det)) return 0;
m_inv[0] = (m[1]*m[2]-m[3]*m[3])/det;
m_inv[1] = (m[0]*m[2]-m[4]*m[4])/det;
m_inv[2] = (m[0]*m[1]-m[5]*m[5])/det;
m_inv[3] = (m[4]*m[5]-m[0]*m[3])/det;
m_inv[4] = (m[3]*m[5]-m[1]*m[4])/det;
m_inv[5] = (m[3]*m[4]-m[2]*m[5])/det;
return 1; }
inline void form_dot(form &dest, form &src) { // m_a.m_b
form m;
memcpy(m, dest, sizeof(form));
dest[0] = m[0]*src[0]+m[4]*src[4]+m[5]*src[5];
dest[1] = m[1]*src[1]+m[3]*src[3]+m[5]*src[5];
dest[2] = m[2]*src[2]+m[3]*src[3]+m[4]*src[4];
dest[3] = m[3]*src[2]+m[1]*src[3]+m[5]*src[4];
dest[4] = m[4]*src[2]+m[5]*src[3]+m[0]*src[4];
dest[5] = m[5]*src[1]+m[4]*src[3]+m[0]*src[5]; }
inline void form_vec_dot(vector &dest, form &m) { // m.a
vector a;
memcpy(a, dest, sizeof(vector));
dest[0] = m[0]*a[0]+m[5]*a[1]+m[4]*a[2];
dest[1] = m[5]*a[0]+m[1]*a[1]+m[3]*a[2];
dest[2] = m[4]*a[0]+m[3]*a[1]+m[2]*a[2]; }
inline void form_to_lform(lform &dest, form &src) {
dest[0] = (long)src[0]; dest[1] = (long)src[1]; dest[2] = (long)src[2];
dest[3] = (long)src[3]; dest[4] = (long)src[4]; dest[5] = (long)src[5]; }
inline void lform_to_form(form &dest, lform &src) {
dest[0] = (double)src[0]; dest[1] = (double)src[1]; dest[2] = (double)src[2];
dest[3] = (double)src[3]; dest[4] = (double)src[4]; dest[5] = (double)src[5];}
// matrix operators
inline double matrix_det(matrix &m) { // |m|
vector axb;
memcpy(&axb, m[0], sizeof(vector));
vec_cross(axb, m[1]);
return vec_dot(axb, m[2]); }
inline void matrix_null(matrix &dest) {
memset(dest, 0, sizeof(dest)); }
inline void matrix_unit(matrix &dest) {
memset(dest, 0, sizeof(dest));
dest[0][0] = dest[1][1] = dest[2][2] = 1.0; }
inline void matrix_scalar_mult(matrix &dest, double f) { // f*m
vec_scalar_mult(dest[0], f);
vec_scalar_mult(dest[1], f);
vec_scalar_mult(dest[2], f); }
inline void matrix_trans(matrix &dest) { // m^t
double f = dest[0][1]; dest[0][1] = dest[1][0]; dest[1][0] = f;
f = dest[0][2]; dest[0][2] = dest[2][0]; dest[2][0] = f;
f = dest[1][2]; dest[1][2] = dest[2][1]; dest[2][1] = f; }
inline int matrix_inv(matrix &dest, matrix &src) { // m^-1
double f = matrix_det(src);
if (fzero(f)) return 0; // singular matrix
memcpy(dest[0], src[1], sizeof(vector));
memcpy(dest[1], src[2], sizeof(vector));
memcpy(dest[2], src[0], sizeof(vector));
vec_cross(dest[0], src[2]);
vec_cross(dest[1], src[0]);
vec_cross(dest[2], src[1]);
matrix_scalar_mult(dest, 1.0/f);
matrix_trans(dest);
return 0; }
inline void matrix_vec_dot(vector &dest, vector &src, matrix &m) { // m.a
dest[0] = m[0][0]*src[0]+m[1][0]*src[1]+m[2][0]*src[2];
dest[1] = m[0][1]*src[0]+m[1][1]*src[1]+m[2][1]*src[2];
dest[2] = m[0][2]*src[0]+m[1][2]*src[1]+m[2][2]*src[2]; }
inline void matrix_to_shape(shape &dest, matrix &src) { // h = m
dest[0] = src[0][0]; dest[1] = src[1][1]; dest[2] = src[2][2];
dest[3] = src[2][1]; dest[4] = src[2][0]; dest[5] = src[1][0]; }
inline void matrix_to_lmatrix(lmatrix &dest, matrix &src) {
vec_to_lvec(dest[0], src[0]);
vec_to_lvec(dest[1], src[1]);
vec_to_lvec(dest[2], src[2]); }
inline void lmatrix_to_matrix(matrix &dest, lmatrix &src) {
lvec_to_vec(dest[0], src[0]);
lvec_to_vec(dest[1], src[1]);
lvec_to_vec(dest[2], src[2]); }
// quaternion operators
inline double quat_dot(quaternion &p, quaternion &q) { // p.q
return p[0]*q[0]+p[1]*q[1]+p[2]*q[2]+p[3]*q[3];
}
inline void quat_norm(quaternion &q) { // q = q/|q|
double f = sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]);
q[0] /= f; q[1] /= f; q[2] /= f; q[3] /= f;
}
inline void quat_conj(quaternion &q) { // q = conj(q)
q[1] = -q[1]; q[2] = -q[2]; q[3] = -q[3];
}
inline void quat_mult(quaternion &dest, quaternion &src) { // dest *= src
quaternion q;
memcpy(q, dest, sizeof(quaternion));
dest[0] = src[0]*q[3]+src[1]*q[2]-src[2]*q[1]+src[3]*q[0];
dest[1] = -src[0]*q[2]+src[1]*q[3]+src[2]*q[0]+src[3]*q[1];
dest[2] = src[0]*q[1]-src[1]*q[0]+src[2]*q[3]+src[3]*q[2];
dest[3] = -src[0]*q[0]-src[1]*q[1]-src[2]*q[2]+src[3]*q[3];
}
inline void quat_div(quaternion &dest, quaternion &src) { // dest /= src
quaternion q;
memcpy(q, dest, sizeof(quaternion));
dest[0] = src[0]*q[3]-src[1]*q[2]+src[2]*q[1]-src[3]*q[0];
dest[1] = -src[0]*q[2]-src[1]*q[3]-src[2]*q[0]-src[3]*q[1];
dest[2] = src[0]*q[1]+src[1]*q[0]-src[2]*q[3]-src[3]*q[2];
dest[3] = -src[0]*q[0]+src[1]*q[1]+src[2]*q[2]-src[3]*q[3];
}
// dest = q*src*conj(q)
inline void quat_vec_rot(vector &dest, vector &src, quaternion &q) {
quaternion aa={q[0]*q[0], q[1]*q[1], q[2]*q[2], q[3]*q[3]};
form ab={q[0]*q[1], q[0]*q[2], q[0]*q[3], q[1]*q[2], q[1]*q[3], q[2]*q[3]};
dest[0] = (aa[0]+aa[1]-aa[2]-aa[3])*src[0]+
((ab[3]-ab[2])*src[1]+(ab[1]+ab[4])*src[2])*2.0;
dest[1] = (aa[0]-aa[1]+aa[2]-aa[3])*src[1]+
((ab[2]+ab[3])*src[0]+(ab[5]-ab[0])*src[2])*2.0;
dest[2] = (aa[0]-aa[1]-aa[2]+aa[3])*src[2]+
((ab[4]-ab[1])*src[0]+(ab[0]+ab[5])*src[1])*2.0;
}
// vector4 operators
inline void vec4_null(vector4 &dest) {
memset(dest, 0, sizeof(vector4));
}
inline double vec4_dot(vector4 &a, vector4 &b) {
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]+a[3]*b[3]; }
// form operators
inline void form4_null(form4 &dest) {
memset(dest, 0, sizeof(form4)); }
inline void form4_unit(form4 &dest) {
memset(dest, 0, sizeof(form4));
dest[0] = dest[1] = dest[2] = dest[3] = 1.0; }
inline double form4_det(form4 &m) {
register double f = m[6]*m[7]-m[5]*m[8];
return m[0]*(
m[1]*(m[2]*m[3]-m[4]*m[4])+
m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])+f*f+
-m[1]*(m[2]*m[6]*m[6]+m[8]*(m[3]*m[8]-2.0*m[4]*m[6]))+
m[9]*(
2.0*(m[2]*m[5]*m[6]+m[3]*m[7]*m[8]-m[4]*(m[6]*m[7]+m[5]*m[8]))+
m[9]*(m[4]*m[4]-m[2]*m[3])); }
inline int form4_inv(form4 &m_inv, form4 &m) {
register double det = form4_det(m);
if (fzero(det)) return 0;
m_inv[0] = (m[1]*(m[2]*m[3]-m[4]*m[4])+
m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])/det;
m_inv[1] = (m[0]*(m[2]*m[3]-m[4]*m[4])+
m[6]*(2.0*m[4]*m[8]-m[2]*m[6])-m[3]*m[8]*m[8])/det;
m_inv[2] = (m[0]*(m[1]*m[3]-m[5]*m[5])+
m[6]*(2.0*m[5]*m[9]-m[1]*m[6])-m[3]*m[9]*m[9])/det;
m_inv[3] = (m[0]*(m[1]*m[2]-m[7]*m[7])+
m[8]*(2.0*m[7]*m[9]-m[1]*m[8])-m[2]*m[9]*m[9])/det;
m_inv[4] = (m[0]*(m[5]*m[7]-m[1]*m[4])+m[1]*m[6]*m[8]+
m[9]*(m[4]*m[9]-m[6]*m[7]-m[5]*m[8]))/det;
m_inv[5] = (m[0]*(m[4]*m[7]-m[2]*m[5])+m[2]*m[6]*m[9]+
m[8]*(m[5]*m[8]-m[6]*m[7]-m[4]*m[9]))/det;
m_inv[6] = (m[1]*(m[4]*m[8]-m[2]*m[6])+m[2]*m[5]*m[9]+
m[7]*(m[6]*m[7]-m[5]*m[8]-m[4]*m[9]))/det;
m_inv[7] = (m[0]*(m[4]*m[5]-m[3]*m[7])+m[3]*m[8]*m[9]+
m[6]*(m[6]*m[7]-m[5]*m[8]-m[4]*m[9]))/det;
m_inv[8] = (m[1]*(m[4]*m[6]-m[3]*m[8])+m[3]*m[7]*m[9]+
m[5]*(m[5]*m[8]-m[6]*m[7]-m[4]*m[9]))/det;
m_inv[9] = (m[2]*(m[5]*m[6]-m[3]*m[9])+m[3]*m[7]*m[8]+
m[4]*(m[4]*m[9]-m[6]*m[7]-m[5]*m[8]))/det;
return 1; }
inline void form4_vec4_dot(vector4 &dest, form4 &m) {
vector4 a;
memcpy(a, dest, sizeof(vector4));
dest[0] = m[0]*a[0]+m[9]*a[1]+m[7]*a[2]+m[6]*a[3];
dest[1] = m[9]*a[0]+m[1]*a[1]+m[6]*a[2]+m[5]*a[3];
dest[2] = m[8]*a[0]+m[7]*a[1]+m[2]*a[2]+m[4]*a[3];
dest[3] = m[6]*a[0]+m[5]*a[1]+m[4]*a[2]+m[3]*a[3]; }
// square regression: y = eqn[0] + eqn[1]*x + eqn[2]*x*x
inline int regress2(vector &eqn, int order, double *x, double *y, int n) {
form xtx = FORM_NULL, xtx_inv;
vector xty = VECTOR_NULL;
double xn, xi, yi;
int i;
vec_null(eqn);
xtx[0] = n;
if ((order = order%2)<0) order = -order; // max: quad regress
if (order<1) xtx[1] = 1.0;
if (order<2) xtx[2] = 1.0;
for (i=0; i<n; ++i) {
xty[0] += (yi = y[i]);
if (order<1) continue;
xty[1] += yi*(xi = xn = x[i]); xtx[5] += xn; xtx[1] += (xn *= xi);
if (order<2) continue;
xty[2] += yi*xn; xtx[3] += (xn *= xi); xtx[2] += xn*xi;
}
if (order>1) xtx[4] = xtx[2];
if (!form_inv(xtx_inv, xtx)) return 0;
memcpy(eqn, xty, sizeof(vector));
form_vec_dot(eqn, xtx_inv);
return 1; }
// cubic regression: y = eqn[0] + eqn[1]*x + eqn[2]*x*x + eqn[3]*x*x*x
inline int regress3(vector4 &eqn, int order, double *x, double *y, int n) {
form4 xtx = FORM4_NULL, xtx_inv;
vector4 xty = VECTOR4_NULL;
double xn, xi, yi;
int i;
vec4_null(eqn);
xtx[0] = n;
if ((order = order%3)<0) order = -order; // max: cubic regress
if (order<1) xtx[1] = 1.0;
if (order<2) xtx[2] = 1.0;
if (order<3) xtx[3] = 1.0;
for (i=0; i<n; ++i) {
xty[0] += (yi = y[i]);
if (order<1) continue;
xty[1] += yi*(xi = xn = x[i]); xtx[9] += xn; xtx[1] += (xn *= xi);
if (order<2) continue;
xty[2] += yi*xn; xtx[7] += (xn *= xi); xtx[2] += xn*xi;
if (order<3) continue;
xty[3] += yi*xn; xtx[4] += (xn *= xi*xi); xtx[3] += xn*xi;
}
if (order>1) xtx[8] = xtx[1];
if (order>2) { xtx[6] = xtx[7]; xtx[5] = xtx[2]; }
if (!form4_inv(xtx_inv, xtx)) return 0;
memcpy(eqn, xty, sizeof(vector4));
form4_vec4_dot(eqn, xtx_inv);
return 1; }
}
#endif

1168
src/USER-EWALDN/pair_buck_coul.cpp
View File

File diff suppressed because it is too large Load Diff

76
src/USER-EWALDN/pair_buck_coul.h
View File

@ -1,76 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(buck/coul,PairBuckCoul)
#else
#ifndef LMP_PAIR_BUCK_COUL_H
#define LMP_PAIR_BUCK_COUL_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBuckCoul : public Pair {
public:
double cut_coul;
PairBuckCoul(class LAMMPS *);
~PairBuckCoul();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
protected:
double cut_buck_global;
double **cut_buck, **cut_buck_read, **cut_bucksq;
double cut_coulsq;
double **buck_a_read, **buck_a, **buck_c_read, **buck_c;
double **buck1, **buck2, **buck_rho_read, **buck_rho, **rhoinv, **offset;
double *cut_respa;
double g_ewald;
int ewald_order, ewald_off;
double tabinnersq;
double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable;
double *etable, *detable, *ptable, *dptable, *vtable, *dvtable;
int ncoulshiftbits, ncoulmask;
void options(char **arg, int order);
void allocate();
void init_tables();
void free_tables();
};
}
#endif
#endif

1160
src/USER-EWALDN/pair_lj_coul.cpp
View File

File diff suppressed because it is too large Load Diff

75
src/USER-EWALDN/pair_lj_coul.h
View File

@ -1,75 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/coul,PairLJCoul)
#else
#ifndef LMP_PAIR_LJ_COUL_H
#define LMP_PAIR_LJ_COUL_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCoul : public Pair {
public:
double cut_coul;
PairLJCoul(class LAMMPS *);
virtual ~PairLJCoul();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
protected:
double cut_lj_global;
double **cut_lj, **cut_lj_read, **cut_ljsq;
double cut_coulsq;
double **epsilon_read, **epsilon, **sigma_read, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset;
double *cut_respa;
double g_ewald;
int ewald_order, ewald_off;
double tabinnersq;
double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable;
double *etable, *detable, *ptable, *dptable, *vtable, *dvtable;
int ncoulshiftbits, ncoulmask;
void options(char **arg, int order);
void allocate();
void init_tables();
void free_tables();
};
}
#endif
#endif

26
src/kspace.cpp
View File

@ -38,6 +38,11 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
order = 5; order = 5;
gridflag = 0; gridflag = 0;
gewaldflag = 0; gewaldflag = 0;
order_6 = 5;
gridflag_6 = 0;
gewaldflag_6 = 0;
slabflag = 0; slabflag = 0;
differentiation_flag = 0; differentiation_flag = 0;
slab_volfactor = 1; slab_volfactor = 1;
@ -272,12 +277,23 @@ void KSpace::modify_params(int narg, char **arg)
if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0) gridflag = 0; if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0) gridflag = 0;
else gridflag = 1; else gridflag = 1;
iarg += 4; iarg += 4;
} else if (strcmp(arg[iarg],"mesh/disp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
nx_pppm_6 = atoi(arg[iarg+1]);
ny_pppm_6 = atoi(arg[iarg+2]);
nz_pppm_6 = atoi(arg[iarg+3]);
if (nx_pppm_6 == 0 || ny_pppm_6 == 0 || nz_pppm_6 == 0) gridflag_6 = 0;
else gridflag_6 = 1;
iarg += 4;
} else if (strcmp(arg[iarg],"order") == 0) { } else if (strcmp(arg[iarg],"order") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order = atoi(arg[iarg+1]); order = atoi(arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"splitorder") == 0) { } else if (strcmp(arg[iarg],"order/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
order_6 = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg],"order/split") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
split_order = atoi(arg[iarg+1]); split_order = atoi(arg[iarg+1]);
if (split_order < 2 || split_order > 6) if (split_order < 2 || split_order > 6)
error->all(FLERR,"Illegal kspace_modify command"); error->all(FLERR,"Illegal kspace_modify command");
@ -292,6 +308,12 @@ void KSpace::modify_params(int narg, char **arg)
if (g_ewald == 0.0) gewaldflag = 0; if (g_ewald == 0.0) gewaldflag = 0;
else gewaldflag = 1; else gewaldflag = 1;
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"gewald/disp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
g_ewald_6 = atof(arg[iarg+1]);
if (g_ewald_6 == 0.0) gewaldflag_6 = 0;
else gewaldflag_6 = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"slab") == 0) { } else if (strcmp(arg[iarg],"slab") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
if (strcmp(arg[iarg+1],"nozforce") == 0) { if (strcmp(arg[iarg+1],"nozforce") == 0) {

11
src/kspace.h
View File

@ -23,14 +23,15 @@ class KSpace : protected Pointers {
friend class FixOMP; friend class FixOMP;
public: public:
double energy; // accumulated energy double energy; // accumulated energy
double energy_1,energy_6;
double virial[6]; // accumlated virial double virial[6]; // accumlated virial
double *eatom,**vatom; // accumulated per-atom energy/virial double *eatom,**vatom; // accumulated per-atom energy/virial
double e2group; // accumulated group-group energy double e2group; // accumulated group-group energy
double f2group[3]; // accumulated group-group force double f2group[3]; // accumulated group-group force
double g_ewald; double g_ewald,g_ewald_6;
int nx_pppm,ny_pppm,nz_pppm; int nx_pppm,ny_pppm,nz_pppm;
int nx_pppm_6,ny_pppm_6,nz_pppm_6;
int nx_msm_max,ny_msm_max,nz_msm_max; int nx_msm_max,ny_msm_max,nz_msm_max;
int group_group_enable; // 1 if style supports group/group calculation int group_group_enable; // 1 if style supports group/group calculation
@ -61,8 +62,10 @@ class KSpace : protected Pointers {
double dgamma(const double &); double dgamma(const double &);
protected: protected:
int gridflag,gewaldflag,differentiation_flag; int gridflag,gridflag_6;
int order,split_order; int gewaldflag,gewaldflag_6;
int order,order_6,split_order;
int differentiation_flag;
int slabflag; int slabflag;
int suffix_flag; // suffix compatibility flag int suffix_flag; // suffix compatibility flag
double scale; double scale;

2
src/pair_coul_dsf.cpp
View File

@ -327,7 +327,7 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void *PairCoulDSF::extract(char *str, int &dim) void *PairCoulDSF::extract(const char *str, int &dim)
{ {
if (strcmp(str,"cut_coul") == 0) { if (strcmp(str,"cut_coul") == 0) {
dim = 0; dim = 0;

2
src/pair_coul_dsf.h
View File

@ -38,7 +38,7 @@ class PairCoulDSF : public Pair {
void write_restart_settings(FILE *); void write_restart_settings(FILE *);
void read_restart_settings(FILE *); void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &); double single(int, int, int, int, double, double, double, double &);
void *extract(char *, int &); void *extract(const char *, int &);
protected: protected:
double cut_coul,cut_coulsq; double cut_coul,cut_coulsq;

2
src/pair_lj_cut_coul_dsf.cpp
View File

@ -460,7 +460,7 @@ double PairLJCutCoulDSF::single(int i, int j, int itype, int jtype, double rsq,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void *PairLJCutCoulDSF::extract(char *str, int &dim) void *PairLJCutCoulDSF::extract(const char *str, int &dim)
{ {
if (strcmp(str,"cut_coul") == 0) { if (strcmp(str,"cut_coul") == 0) {
dim = 0; dim = 0;

2
src/pair_lj_cut_coul_dsf.h
View File

@ -38,7 +38,7 @@ class PairLJCutCoulDSF : public Pair {
void write_restart_settings(FILE *); void write_restart_settings(FILE *);
void read_restart_settings(FILE *); void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &); double single(int, int, int, int, double, double, double, double &);
void *extract(char *, int &); void *extract(const char *, int &);
protected: protected:
double cut_lj_global; double cut_lj_global;