Skip to content

Commit

Permalink
Merge pull request #46 from rk-lindsey/develop
Browse files Browse the repository at this point in the history 
Merge in BB optimize2 and fix NPT bug for LAMMPS linking
  • Loading branch information
@rk-lindsey
rk-lindsey authored Jan 5, 2023
2 parents 1a628cc + 1a547e0 commit ed5a315
Show file tree
Hide file tree
Showing 32 changed files with 12,492 additions and 2,616 deletions.
276 changes: 126 additions & 150 deletions chimesFF/api/chimescalc_C.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,9 @@ void chimes_compute_2b_props(double rij, double dr[3], char *atype2b[2], double
// convert all doubles, etc., from C to type vector for C++
// declare needed vectors for chimes
vector <double> dr_vec(3);

chimes2BTmp chimes_tmp(chimes_ptr->poly_orders[0]) ;

dr_vec[0] = dr[0];
dr_vec[1] = dr[1];
dr_vec[2] = dr[2];
Expand All @@ -124,41 +127,29 @@ void chimes_compute_2b_props(double rij, double dr[3], char *atype2b[2], double
}
//type_vec[0] = chimes_ptr->atmtoidx[atype2b[0]];
//type_vec[1] = chimes_ptr->atmtoidx[atype2b[1]];
vector<vector<double*> >force_vec;
force_vec.resize(2,vector<double*>(3));
force_vec[0][0] = &force[0][0];
force_vec[0][1] = &force[0][1];
force_vec[0][2] = &force[0][2];
force_vec[1][0] = &force[1][0];
force_vec[1][1] = &force[1][1];
force_vec[1][2] = &force[1][2];
vector<double*> stress_vec(9);
stress_vec[0] = &stress[0];
stress_vec[1] = &stress[1];
stress_vec[2] = &stress[2];
stress_vec[3] = &stress[3];
stress_vec[4] = &stress[4];
stress_vec[5] = &stress[5];
stress_vec[6] = &stress[6];
stress_vec[7] = &stress[7];
stress_vec[8] = &stress[8];
chimes_ptr->compute_2B(rij, dr_vec, type_vec, force_vec, stress_vec, *epot);
vector<double> force_vec(2*CHDIM,0.0) ;;
vector<double> stress_vec(9,0.0);

chimes_ptr->compute_2B(rij, dr_vec, type_vec, force_vec, stress_vec, *epot, chimes_tmp);


// save forces and stress tensor components
force[0][0] = *force_vec[0][0];
force[0][1] = *force_vec[0][1];
force[0][2] = *force_vec[0][2];
force[1][0] = *force_vec[1][0];
force[1][1] = *force_vec[1][1];
force[1][2] = *force_vec[1][2];
stress[0] = *stress_vec[0];
stress[1] = *stress_vec[1];
stress[2] = *stress_vec[2];
stress[3] = *stress_vec[3];
stress[4] = *stress_vec[4];
stress[5] = *stress_vec[5];
stress[6] = *stress_vec[6];
stress[7] = *stress_vec[7];
stress[8] = *stress_vec[8];
force[0][0] += force_vec[0*CHDIM+0];
force[0][1] += force_vec[0*CHDIM+1];
force[0][2] += force_vec[0*CHDIM+2];
force[1][0] += force_vec[1*CHDIM+0];
force[1][1] += force_vec[1*CHDIM+1];
force[1][2] += force_vec[1*CHDIM+2];

stress[0] += stress_vec[0];
stress[1] += stress_vec[1];
stress[2] += stress_vec[2];
stress[3] += stress_vec[3];
stress[4] += stress_vec[4];
stress[5] += stress_vec[5];
stress[6] += stress_vec[6];
stress[7] += stress_vec[7];
stress[8] += stress_vec[8];
}

void chimes_compute_3b_props(double dr_3b[3], double dist_3b[3][3], char *atype3b[3], double f3b[3][3], double stress[9], double *epot) {
Expand All @@ -167,20 +158,24 @@ void chimes_compute_3b_props(double dr_3b[3], double dist_3b[3][3], char *atype3
dr_3b_vec[0] = dr_3b[0];
dr_3b_vec[1] = dr_3b[1];
dr_3b_vec[2] = dr_3b[2];
vector< vector<double> > dist_3b_vec(3, vector<double> (3));
dist_3b_vec[0][0] = dist_3b[0][0];
dist_3b_vec[0][1] = dist_3b[0][1];
dist_3b_vec[0][2] = dist_3b[0][2];
dist_3b_vec[1][0] = dist_3b[1][0];
dist_3b_vec[1][1] = dist_3b[1][1];
dist_3b_vec[1][2] = dist_3b[1][2];
dist_3b_vec[2][0] = dist_3b[2][0];
dist_3b_vec[2][1] = dist_3b[2][1];
dist_3b_vec[2][2] = dist_3b[2][2];

vector<double> dist_3b_vec(CHDIM*3) ;

dist_3b_vec[0*CHDIM+0] = dist_3b[0][0];
dist_3b_vec[0*CHDIM+1] = dist_3b[0][1];
dist_3b_vec[0*CHDIM+2] = dist_3b[0][2];
dist_3b_vec[1*CHDIM+0] = dist_3b[1][0];
dist_3b_vec[1*CHDIM+1] = dist_3b[1][1];
dist_3b_vec[1*CHDIM+2] = dist_3b[1][2];
dist_3b_vec[2*CHDIM+0] = dist_3b[2][0];
dist_3b_vec[2*CHDIM+1] = dist_3b[2][1];
dist_3b_vec[2*CHDIM+2] = dist_3b[2][2];

vector <int> type_3b_vec(3);
type_3b_vec[0] = distance(chimes_ptr->atmtyps.begin(),find(chimes_ptr->atmtyps.begin(), chimes_ptr->atmtyps.end(), atype3b[0]));
type_3b_vec[1] = distance(chimes_ptr->atmtyps.begin(),find(chimes_ptr->atmtyps.begin(), chimes_ptr->atmtyps.end(), atype3b[1]));
type_3b_vec[2] = distance(chimes_ptr->atmtyps.begin(),find(chimes_ptr->atmtyps.begin(), chimes_ptr->atmtyps.end(), atype3b[2]));

for(int i=0; i<3; i++)
{
if (type_3b_vec[i] >= chimes_ptr->atmtyps.size())
Expand All @@ -192,50 +187,38 @@ void chimes_compute_3b_props(double dr_3b[3], double dist_3b[3][3], char *atype3
//type_3b_vec[0] = chimes_ptr->atmtoidx[atype3b[0]];
//type_3b_vec[1] = chimes_ptr->atmtoidx[atype3b[1]];
//type_3b_vec[2] = chimes_ptr->atmtoidx[atype3b[2]];
vector<vector<double*> >force_3b_vec;
force_3b_vec.resize(3,vector<double*>(3));
force_3b_vec[0][0] = &f3b[0][0];
force_3b_vec[0][1] = &f3b[0][1];
force_3b_vec[0][2] = &f3b[0][2];
force_3b_vec[1][0] = &f3b[1][0];
force_3b_vec[1][1] = &f3b[1][1];
force_3b_vec[1][2] = &f3b[1][2];
force_3b_vec[2][0] = &f3b[2][0];
force_3b_vec[2][1] = &f3b[2][1];
force_3b_vec[2][2] = &f3b[2][2];
vector<double*> stress_vec(9);
stress_vec[0] = &stress[0];
stress_vec[1] = &stress[1];
stress_vec[2] = &stress[2];
stress_vec[3] = &stress[3];
stress_vec[4] = &stress[4];
stress_vec[5] = &stress[5];
stress_vec[6] = &stress[6];
stress_vec[7] = &stress[7];
stress_vec[8] = &stress[8];
chimes_ptr->compute_3B(dr_3b_vec, dist_3b_vec, type_3b_vec, force_3b_vec, stress_vec, *epot);
vector<double> force_3b_vec(3*CHDIM, 0.0) ;
vector<double> stress_vec(9,0.0);
chimes3BTmp chimes_tmp(chimes_ptr->poly_orders[1]) ;

chimes_ptr->compute_3B(dr_3b_vec, dist_3b_vec, type_3b_vec, force_3b_vec, stress_vec, *epot,
chimes_tmp) ;

// save forces and stress tensor components
f3b[0][0] = *force_3b_vec[0][0];
f3b[0][1] = *force_3b_vec[0][1];
f3b[0][2] = *force_3b_vec[0][2];
f3b[1][0] = *force_3b_vec[1][0];
f3b[1][1] = *force_3b_vec[1][1];
f3b[1][2] = *force_3b_vec[1][2];
f3b[2][0] = *force_3b_vec[2][0];
f3b[2][1] = *force_3b_vec[2][1];
f3b[2][2] = *force_3b_vec[2][2];
stress[0] = *stress_vec[0];
stress[1] = *stress_vec[1];
stress[2] = *stress_vec[2];
stress[3] = *stress_vec[3];
stress[4] = *stress_vec[4];
stress[5] = *stress_vec[5];
stress[6] = *stress_vec[6];
stress[7] = *stress_vec[7];
stress[8] = *stress_vec[8];

f3b[0][0] += force_3b_vec[0*CHDIM+0];
f3b[0][1] += force_3b_vec[0*CHDIM+1];
f3b[0][2] += force_3b_vec[0*CHDIM+2];
f3b[1][0] += force_3b_vec[1*CHDIM+0];
f3b[1][1] += force_3b_vec[1*CHDIM+1];
f3b[1][2] += force_3b_vec[1*CHDIM+2];
f3b[2][0] += force_3b_vec[2*CHDIM+0];
f3b[2][1] += force_3b_vec[2*CHDIM+1];
f3b[2][2] += force_3b_vec[2*CHDIM+2];

stress[0] += stress_vec[0];
stress[1] += stress_vec[1];
stress[2] += stress_vec[2];
stress[3] += stress_vec[3];
stress[4] += stress_vec[4];
stress[5] += stress_vec[5];
stress[6] += stress_vec[6];
stress[7] += stress_vec[7];
stress[8] += stress_vec[8];
}

void chimes_compute_4b_props(double dr_4b[6], double dist_4b[6][3], char *atype4b[4], double f4b[4][3], double stress[9], double *epot) {
void chimes_compute_4b_props(double dr_4b[6], double dist_4b[6][3], char *atype4b[4], double f4b[4][3], double stress[9], double *epot)
{
// convert all doubles, etc., from C to type vector for C++
vector <double> dr_4b_vec(6);
dr_4b_vec[0] = dr_4b[0];
Expand All @@ -244,25 +227,28 @@ void chimes_compute_4b_props(double dr_4b[6], double dist_4b[6][3], char *atype4
dr_4b_vec[3] = dr_4b[3];
dr_4b_vec[4] = dr_4b[4];
dr_4b_vec[5] = dr_4b[5];
vector< vector<double> > dist_4b_vec(6, vector<double> (3));
dist_4b_vec[0][0] = dist_4b[0][0];
dist_4b_vec[0][1] = dist_4b[0][1];
dist_4b_vec[0][2] = dist_4b[0][2];
dist_4b_vec[1][0] = dist_4b[1][0];
dist_4b_vec[1][1] = dist_4b[1][1];
dist_4b_vec[1][2] = dist_4b[1][2];
dist_4b_vec[2][0] = dist_4b[2][0];
dist_4b_vec[2][1] = dist_4b[2][1];
dist_4b_vec[2][2] = dist_4b[2][2];
dist_4b_vec[3][0] = dist_4b[3][0];
dist_4b_vec[3][1] = dist_4b[3][1];
dist_4b_vec[3][2] = dist_4b[3][2];
dist_4b_vec[4][0] = dist_4b[4][0];
dist_4b_vec[4][1] = dist_4b[4][1];
dist_4b_vec[4][2] = dist_4b[4][2];
dist_4b_vec[5][0] = dist_4b[5][0];
dist_4b_vec[5][1] = dist_4b[5][1];
dist_4b_vec[5][2] = dist_4b[5][2];

vector<double> dist_4b_vec(6*CHDIM);

dist_4b_vec[0*CHDIM+0] = dist_4b[0][0];
dist_4b_vec[0*CHDIM+1] = dist_4b[0][1];
dist_4b_vec[0*CHDIM+2] = dist_4b[0][2];
dist_4b_vec[1*CHDIM+0] = dist_4b[1][0];
dist_4b_vec[1*CHDIM+1] = dist_4b[1][1];
dist_4b_vec[1*CHDIM+2] = dist_4b[1][2];
dist_4b_vec[2*CHDIM+0] = dist_4b[2][0];
dist_4b_vec[2*CHDIM+1] = dist_4b[2][1];
dist_4b_vec[2*CHDIM+2] = dist_4b[2][2];
dist_4b_vec[3*CHDIM+0] = dist_4b[3][0];
dist_4b_vec[3*CHDIM+1] = dist_4b[3][1];
dist_4b_vec[3*CHDIM+2] = dist_4b[3][2];
dist_4b_vec[4*CHDIM+0] = dist_4b[4][0];
dist_4b_vec[4*CHDIM+1] = dist_4b[4][1];
dist_4b_vec[4*CHDIM+2] = dist_4b[4][2];
dist_4b_vec[5*CHDIM+0] = dist_4b[5][0];
dist_4b_vec[5*CHDIM+1] = dist_4b[5][1];
dist_4b_vec[5*CHDIM+2] = dist_4b[5][2];

vector <int> type_4b_vec(4);
type_4b_vec[0] = distance(chimes_ptr->atmtyps.begin(),find(chimes_ptr->atmtyps.begin(), chimes_ptr->atmtyps.end(), atype4b[0]));
type_4b_vec[1] = distance(chimes_ptr->atmtyps.begin(),find(chimes_ptr->atmtyps.begin(), chimes_ptr->atmtyps.end(), atype4b[1]));
Expand All @@ -280,51 +266,41 @@ void chimes_compute_4b_props(double dr_4b[6], double dist_4b[6][3], char *atype4
//type_4b_vec[1] = chimes_ptr->atmtoidx[atype4b[1]];
//type_4b_vec[2] = chimes_ptr->atmtoidx[atype4b[2]];
//type_4b_vec[3] = chimes_ptr->atmtoidx[atype4b[3]];
vector<vector<double*> >force_4b_vec;
force_4b_vec.resize(4,vector<double*>(3));
force_4b_vec[0][0] = &f4b[0][0];
force_4b_vec[0][1] = &f4b[0][1];
force_4b_vec[0][2] = &f4b[0][2];
force_4b_vec[1][0] = &f4b[1][0];
force_4b_vec[1][1] = &f4b[1][1];
force_4b_vec[1][2] = &f4b[1][2];
force_4b_vec[2][0] = &f4b[2][0];
force_4b_vec[2][1] = &f4b[2][1];
force_4b_vec[2][2] = &f4b[2][2];
force_4b_vec[3][0] = &f4b[3][0];
force_4b_vec[3][1] = &f4b[3][1];
force_4b_vec[3][2] = &f4b[3][2];
vector<double*> stress_vec(9);
stress_vec[0] = &stress[0];
stress_vec[1] = &stress[1];
stress_vec[2] = &stress[2];
stress_vec[3] = &stress[3];
stress_vec[4] = &stress[4];
stress_vec[5] = &stress[5];
stress_vec[6] = &stress[6];
stress_vec[7] = &stress[7];
stress_vec[8] = &stress[8];
chimes_ptr->compute_4B(dr_4b_vec, dist_4b_vec, type_4b_vec, force_4b_vec, stress_vec, *epot);

vector<double> force_4b(4*CHDIM,0.0) ;

// No auto-conversion from double[9] to vector.
vector<double> stress_vec(9) ;

for ( int j = 0 ; j < 9 ; j++ ) {
stress_vec[j] = stress[j] ;
}

chimes4BTmp chimes_tmp(chimes_ptr->poly_orders[2]) ;

chimes_ptr->compute_4B(dr_4b_vec, dist_4b_vec, type_4b_vec, force_4b, stress_vec, *epot, chimes_tmp);
// save forces and stress tensor components
f4b[0][0] = *force_4b_vec[0][0];
f4b[0][1] = *force_4b_vec[0][1];
f4b[0][2] = *force_4b_vec[0][2];
f4b[1][0] = *force_4b_vec[1][0];
f4b[1][1] = *force_4b_vec[1][1];
f4b[1][2] = *force_4b_vec[1][2];
f4b[2][0] = *force_4b_vec[2][0];
f4b[2][1] = *force_4b_vec[2][1];
f4b[2][2] = *force_4b_vec[2][2];
f4b[3][0] = *force_4b_vec[3][0];
f4b[3][1] = *force_4b_vec[3][1];
f4b[3][2] = *force_4b_vec[3][2];
stress[0] = *stress_vec[0];
stress[1] = *stress_vec[1];
stress[2] = *stress_vec[2];
stress[3] = *stress_vec[3];
stress[4] = *stress_vec[4];
stress[5] = *stress_vec[5];
stress[6] = *stress_vec[6];
stress[7] = *stress_vec[7];
stress[8] = *stress_vec[8];

f4b[0][0] += force_4b[0*CHDIM+0];
f4b[0][1] += force_4b[0*CHDIM+1];
f4b[0][2] += force_4b[0*CHDIM+2];
f4b[1][0] += force_4b[1*CHDIM+0];
f4b[1][1] += force_4b[1*CHDIM+1];
f4b[1][2] += force_4b[1*CHDIM+2];
f4b[2][0] += force_4b[2*CHDIM+0];
f4b[2][1] += force_4b[2*CHDIM+1];
f4b[2][2] += force_4b[2*CHDIM+2];
f4b[3][0] += force_4b[3*CHDIM+0];
f4b[3][1] += force_4b[3*CHDIM+1];
f4b[3][2] += force_4b[3*CHDIM+2];

stress[0] = stress_vec[0];
stress[1] = stress_vec[1];
stress[2] = stress_vec[2];
stress[3] = stress_vec[3];
stress[4] = stress_vec[4];
stress[5] = stress_vec[5];
stress[6] = stress_vec[6];
stress[7] = stress_vec[7];
stress[8] = stress_vec[8];
}
24 changes: 12 additions & 12 deletions chimesFF/examples/c/Makefile
View file
Original file line number Diff line number Diff line change
@@ -1,34 +1,34 @@
CC_LOC = $(realpath .)

CXX=g++ -O3 -std=c++11
CXX=g++ -O3 -std=c++11

CHIMESFF_SRC=$(CC_LOC)/../../src/chimesFF.cpp
CHIMESFF_HDR=$(CC_LOC)/../../src/chimesFF.h

chimesFF.o : $(CHIMESFF_SRC)
$(CXX) -c $(CHIMESFF_SRC)

WRAPPER_SRC=$(CC_LOC)/../../api/chimescalc_C.cpp
WRAPPER_HDR=$(CC_LOC)/../../api/chimescalc_C.h

chimescalc_C.o : $(WRAPPER_SRC) $(WRAPPER_HDR) $(CHIMESFF_SRC) $(CHIMESFF_HDR)
$(CXX) -c $(WRAPPER_SRC) $(CHIMESFF_SRC) -I $(CC_LOC)/../../api/ -I $(CC_LOC)/../../src/

CC = gcc

test_wrapper-C.o : main.c $(WRAPPER_SRC) $(WRAPPER_HDR)
$(CC) -c main.c -o test_wrapper-C.o -I $(CC_LOC)/../../api/
$(CC) -c main.c -o test_wrapper-C.o -I $(CC_LOC)/../../api/

LINKS = chimesFF.o test_wrapper-C.o chimescalc_C.o

test_wrapper-C : $(LINKS)
$(CXX) $(LINKS) -o chimescalc-test_direct-C

test_wrapper-C : $(LINKS)
$(CXX) $(LINKS) -o C_wrapper-direct_interface
clean:
rm -f *.o

rm -f *.o
clean-all:
rm -f *.o
rm -f *.o
rm -f chimescalc-test_direct-C

all:
Expand Down
Loading

0 comments on commit ed5a315

Please sign in to comment.