Skip to content

Commit f5835ef

Browse files
amcadmusamcadmus
authored
Merge pull request #142 from amcadmus/devel
support nopbc, fix bug compiling with float_prec=low
2 parents 74f5ece + b88c832 commit f5835ef

File tree

13 files changed

+465
-187
lines changed

13 files changed

+465
-187
lines changed

‎source/lib/include/Ewald.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@ rec_err_esti(const VALUETYPE & test_q,
8181
template <typename VALUETYPE>
8282
void
8383
cmpt_k(vector<int> & KK,
84-
const SimulationRegion<VALUETYPE>& region,
84+
const SimulationRegion<double>& region,
8585
const EwaldParameters<VALUETYPE>& param)
8686
{
8787
const double * boxt_ = region.getBoxTensor();
@@ -154,7 +154,8 @@ EwaldReciprocal(VALUETYPE & ener,
154154
for (int ii = 0; ii < natoms; ++ii){
155155
int thread_id = omp_get_thread_num();
156156
double ir[3];
157-
region.phys2Inter(ir, &coord[ii*3]);
157+
double tmpcoord[3] = {coord[ii*3], coord[ii*3+1], coord[ii*3+2]};
158+
region.phys2Inter(ir, tmpcoord);
158159
for (int mm0 = -KK[0]/2; mm0 <= KK[0]/2; ++mm0){
159160
double mr[3];
160161
mr[0] = ir[0] * mm0;

‎source/lib/include/NeighborList.h‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#include "MathUtilities.h"
88
#include "SimulationRegion.h"
99

10+
// build nlist by an extended grid
1011
void
1112
build_nlist (vector<vector<int > > & nlist0,
1213
vector<vector<int > > & nlist1,
@@ -20,6 +21,8 @@ build_nlist (vector<vector<int > > & nlist0,
2021
const vector<int > & ext_end_,
2122
const SimulationRegion<double> & region,
2223
const vector<int > & global_grid);
24+
25+
// build nlist by a grid for a periodic region
2326
void
2427
build_nlist (vector<vector<int > > & nlist0,
2528
vector<vector<int > > & nlist1,
@@ -28,6 +31,8 @@ build_nlist (vector<vector<int > > & nlist0,
2831
const double & rc1,
2932
const vector<int > & grid,
3033
const SimulationRegion<double> & region);
34+
35+
// build nlist by a grid for a periodic region, atoms selected by sel0 and sel1
3136
void
3237
build_nlist (vector<vector<int > > & nlist0,
3338
vector<vector<int > > & nlist1,
@@ -38,6 +43,19 @@ build_nlist (vector<vector<int > > & nlist0,
3843
const double & rc1,
3944
const vector<int > & grid,
4045
const SimulationRegion<double> & region);
46+
47+
// brute force (all-to-all distance computation) neighbor list building
48+
// if region is NULL, open boundary is assumed,
49+
// otherwise, periodic boundary condition is defined by region
50+
void
51+
build_nlist (vector<vector<int > > & nlist0,
52+
vector<vector<int > > & nlist1,
53+
const vector<double > & coord,
54+
const double & rc0_,
55+
const double & rc1_,
56+
const SimulationRegion<double > * region = NULL);
57+
58+
// copy periodic images for the system
4159
void
4260
copy_coord (vector<double > & out_c,
4361
vector<int > & out_t,

‎source/lib/src/NeighborList.cpp‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -572,6 +572,58 @@ build_nlist (vector<vector<int > > & nlist0,
572572
}
573573
}
574574

575+
576+
void
577+
build_nlist (vector<vector<int > > & nlist0,
578+
vector<vector<int > > & nlist1,
579+
const vector<double > & posi3,
580+
const double & rc0_,
581+
const double & rc1_,
582+
const SimulationRegion<double > * region)
583+
{
584+
double rc0 (rc0_);
585+
double rc1 (rc1_);
586+
assert (rc0 <= rc1);
587+
double rc02 = rc0 * rc0;
588+
// negative rc0 means not applying rc0
589+
if (rc0 < 0) rc02 = 0;
590+
double rc12 = rc1 * rc1;
591+
592+
unsigned natoms = posi3.size()/3;
593+
nlist0.clear();
594+
nlist1.clear();
595+
nlist0.resize(natoms);
596+
nlist1.resize(natoms);
597+
for (unsigned ii = 0; ii < natoms; ++ii){
598+
nlist0[ii].reserve (60);
599+
nlist1[ii].reserve (60);
600+
}
601+
for (unsigned ii = 0; ii < natoms; ++ii){
602+
for (unsigned jj = ii+1; jj < natoms; ++jj){
603+
double diff[3];
604+
if (region != NULL) {
605+
region->diffNearestNeighbor (posi3[jj*3+0], posi3[jj*3+1], posi3[jj*3+2],
606+
posi3[ii*3+0], posi3[ii*3+1], posi3[ii*3+2],
607+
diff[0], diff[1], diff[2]);
608+
}
609+
else {
610+
diff[0] = posi3[jj*3+0] - posi3[ii*3+0];
611+
diff[1] = posi3[jj*3+1] - posi3[ii*3+1];
612+
diff[2] = posi3[jj*3+2] - posi3[ii*3+2];
613+
}
614+
double r2 = MathUtilities::dot<double> (diff, diff);
615+
if (r2 < rc02) {
616+
nlist0[ii].push_back (jj);
617+
nlist0[jj].push_back (ii);
618+
}
619+
else if (r2 < rc12) {
620+
nlist1[ii].push_back (jj);
621+
nlist1[jj].push_back (ii);
622+
}
623+
}
624+
}
625+
}
626+
575627
static int compute_pbc_shift (int idx,
576628
int ncell)
577629
{

‎source/op/descrpt.cc‎

Lines changed: 24 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -120,15 +120,34 @@ class DescrptOp : public OpKernel {
120120

121121
int nei_mode = 0;
122122
if (mesh_tensor.shape().dim_size(0) == 16) {
123+
// lammps neighbor list
123124
nei_mode = 3;
124125
}
125126
else if (mesh_tensor.shape().dim_size(0) == 12) {
127+
// user provided extended mesh
126128
nei_mode = 2;
127129
}
128130
else if (mesh_tensor.shape().dim_size(0) == 6) {
131+
// manual copied pbc
129132
assert (nloc == nall);
130133
nei_mode = 1;
131134
}
135+
else if (mesh_tensor.shape().dim_size(0) == 0) {
136+
// no pbc
137+
nei_mode = -1;
138+
}
139+
else {
140+
throw runtime_error("invalid mesh tensor");
141+
}
142+
bool b_pbc = true;
143+
// if region is given extended, do not use pbc
144+
if (nei_mode >= 1 || nei_mode == -1) {
145+
b_pbc = false;
146+
}
147+
bool b_norm_atom = false;
148+
if (nei_mode == 1){
149+
b_norm_atom = true;
150+
}
132151

133152
// Create an output tensor
134153
TensorShape descrpt_shape ;
@@ -200,7 +219,7 @@ class DescrptOp : public OpKernel {
200219
for (int dd = 0; dd < 3; ++dd){
201220
d_coord3[ii*3+dd] = coord(kk, ii*3+dd);
202221
}
203-
if (nei_mode <= 1){
222+
if (b_norm_atom){
204223
compute_t inter[3];
205224
region.phys2Inter (inter, &d_coord3[3*ii]);
206225
for (int dd = 0; dd < 3; ++dd){
@@ -263,14 +282,11 @@ class DescrptOp : public OpKernel {
263282
}
264283
::build_nlist (d_nlist_a, d_nlist_r, d_coord3, nloc, rcut_a, rcut_r, nat_stt, ncell, ext_stt, ext_end, region, ncell);
265284
}
266-
else {
267-
build_nlist (d_nlist_a, d_nlist_r, rcut_a, rcut_r, d_coord3, region);
285+
else if (nei_mode == -1){
286+
::build_nlist (d_nlist_a, d_nlist_r, d_coord3, rcut_a, rcut_r, NULL);
268287
}
269-
270-
bool b_pbc = true;
271-
// if region is given extended, do not use pbc
272-
if (nei_mode >= 1) {
273-
b_pbc = false;
288+
else {
289+
throw runtime_error("unknow neighbor mode");
274290
}
275291

276292
// loop over atoms, compute descriptors for each atom
@@ -399,48 +415,6 @@ class DescrptOp : public OpKernel {
399415
}
400416
}
401417
void
402-
build_nlist (vector<vector<int > > & nlist0,
403-
vector<vector<int > > & nlist1,
404-
const compute_t & rc0_,
405-
const compute_t & rc1_,
406-
const vector<compute_t > & posi3,
407-
const SimulationRegion<compute_t > & region) const {
408-
compute_t rc0 (rc0_);
409-
compute_t rc1 (rc1_);
410-
assert (rc0 <= rc1);
411-
compute_t rc02 = rc0 * rc0;
412-
// negative rc0 means not applying rc0
413-
if (rc0 < 0) rc02 = 0;
414-
compute_t rc12 = rc1 * rc1;
415-
416-
unsigned natoms = posi3.size()/3;
417-
nlist0.clear();
418-
nlist1.clear();
419-
nlist0.resize(natoms);
420-
nlist1.resize(natoms);
421-
for (unsigned ii = 0; ii < natoms; ++ii){
422-
nlist0[ii].reserve (60);
423-
nlist1[ii].reserve (60);
424-
}
425-
for (unsigned ii = 0; ii < natoms; ++ii){
426-
for (unsigned jj = ii+1; jj < natoms; ++jj){
427-
compute_t diff[3];
428-
region.diffNearestNeighbor (posi3[jj*3+0], posi3[jj*3+1], posi3[jj*3+2],
429-
posi3[ii*3+0], posi3[ii*3+1], posi3[ii*3+2],
430-
diff[0], diff[1], diff[2]);
431-
compute_t r2 = MathUtilities::dot<compute_t> (diff, diff);
432-
if (r2 < rc02) {
433-
nlist0[ii].push_back (jj);
434-
nlist0[jj].push_back (ii);
435-
}
436-
else if (r2 < rc12) {
437-
nlist1[ii].push_back (jj);
438-
nlist1[jj].push_back (ii);
439-
}
440-
}
441-
}
442-
}
443-
void
444418
make_axis (vector<int > & axis_type,
445419
vector<int > & axis_idx,
446420
const int & type,

‎source/op/descrpt_se_a.cc‎

Lines changed: 24 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -119,15 +119,34 @@ class DescrptSeAOp : public OpKernel {
119119

120120
int nei_mode = 0;
121121
if (mesh_tensor.shape().dim_size(0) == 16) {
122+
// lammps neighbor list
122123
nei_mode = 3;
123124
}
124125
else if (mesh_tensor.shape().dim_size(0) == 12) {
126+
// user provided extended mesh
125127
nei_mode = 2;
126128
}
127129
else if (mesh_tensor.shape().dim_size(0) == 6) {
130+
// manual copied pbc
128131
assert (nloc == nall);
129132
nei_mode = 1;
130133
}
134+
else if (mesh_tensor.shape().dim_size(0) == 0) {
135+
// no pbc
136+
nei_mode = -1;
137+
}
138+
else {
139+
throw runtime_error("invalid mesh tensor");
140+
}
141+
bool b_pbc = true;
142+
// if region is given extended, do not use pbc
143+
if (nei_mode >= 1 || nei_mode == -1) {
144+
b_pbc = false;
145+
}
146+
bool b_norm_atom = false;
147+
if (nei_mode == 1){
148+
b_norm_atom = true;
149+
}
131150

132151
// Create an output tensor
133152
TensorShape descrpt_shape ;
@@ -196,7 +215,7 @@ class DescrptSeAOp : public OpKernel {
196215
for (int dd = 0; dd < 3; ++dd){
197216
d_coord3[ii*3+dd] = coord(kk, ii*3+dd);
198217
}
199-
if (nei_mode <= 1){
218+
if (b_norm_atom){
200219
compute_t inter[3];
201220
region.phys2Inter (inter, &d_coord3[3*ii]);
202221
for (int dd = 0; dd < 3; ++dd){
@@ -259,14 +278,11 @@ class DescrptSeAOp : public OpKernel {
259278
}
260279
::build_nlist (d_nlist_a, d_nlist_r, d_coord3, nloc, rcut_a, rcut_r, nat_stt, ncell, ext_stt, ext_end, region, ncell);
261280
}
262-
else {
263-
build_nlist (d_nlist_a, d_nlist_r, rcut_a, rcut_r, d_coord3, region);
281+
else if (nei_mode == -1){
282+
::build_nlist (d_nlist_a, d_nlist_r, d_coord3, rcut_a, rcut_r, NULL);
264283
}
265-
266-
bool b_pbc = true;
267-
// if region is given extended, do not use pbc
268-
if (nei_mode >= 1) {
269-
b_pbc = false;
284+
else {
285+
throw runtime_error("unknow neighbor mode");
270286
}
271287

272288
// loop over atoms, compute descriptors for each atom
@@ -351,48 +367,6 @@ class DescrptSeAOp : public OpKernel {
351367
sec[ii] = sec[ii-1] + n_sel[ii-1];
352368
}
353369
}
354-
void
355-
build_nlist (vector<vector<int > > & nlist0,
356-
vector<vector<int > > & nlist1,
357-
const compute_t & rc0_,
358-
const compute_t & rc1_,
359-
const vector<compute_t > & posi3,
360-
const SimulationRegion<compute_t > & region) const {
361-
compute_t rc0 (rc0_);
362-
compute_t rc1 (rc1_);
363-
assert (rc0 <= rc1);
364-
compute_t rc02 = rc0 * rc0;
365-
// negative rc0 means not applying rc0
366-
if (rc0 < 0) rc02 = 0;
367-
compute_t rc12 = rc1 * rc1;
368-
369-
unsigned natoms = posi3.size()/3;
370-
nlist0.clear();
371-
nlist1.clear();
372-
nlist0.resize(natoms);
373-
nlist1.resize(natoms);
374-
for (unsigned ii = 0; ii < natoms; ++ii){
375-
nlist0[ii].reserve (60);
376-
nlist1[ii].reserve (60);
377-
}
378-
for (unsigned ii = 0; ii < natoms; ++ii){
379-
for (unsigned jj = ii+1; jj < natoms; ++jj){
380-
compute_t diff[3];
381-
region.diffNearestNeighbor (posi3[jj*3+0], posi3[jj*3+1], posi3[jj*3+2],
382-
posi3[ii*3+0], posi3[ii*3+1], posi3[ii*3+2],
383-
diff[0], diff[1], diff[2]);
384-
compute_t r2 = MathUtilities::dot<compute_t> (diff, diff);
385-
if (r2 < rc02) {
386-
nlist0[ii].push_back (jj);
387-
nlist0[jj].push_back (ii);
388-
}
389-
else if (r2 < rc12) {
390-
nlist1[ii].push_back (jj);
391-
nlist1[jj].push_back (ii);
392-
}
393-
}
394-
}
395-
}
396370
};
397371

398372
REGISTER_KERNEL_BUILDER(Name("DescrptSeA").Device(DEVICE_CPU), DescrptSeAOp);

0 commit comments

Comments
 (0)