Skip to content

Commit

Permalink
Merge pull request #56 from LindseyLab-umich/fix-merge-conflict
Browse files Browse the repository at this point in the history
Fix merge conflict
  • Loading branch information
RKLindsey authored Mar 19, 2024
2 parents 231d01d + b923a9a commit 6cb12c0
Show file tree
Hide file tree
Showing 19 changed files with 4,568 additions and 74 deletions.
2 changes: 2 additions & 0 deletions CODEOWNERS
Validating CODEOWNERS rules …
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
* @RKLindsey
* @rk-lindsey
14 changes: 13 additions & 1 deletion chimesFF/api/chimescalc_C.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,16 @@ void init_chimes (int rank) {

void chimes_read_params(char *param_file) {
chimes_ptr->read_parameters(param_file);
chimes_ptr->build_pair_int_trip_map();
chimes_ptr->build_pair_int_quad_map();
}

void chimes_build_pair_int_trip_map() {
chimes_ptr->build_pair_int_trip_map();
}

void chimes_build_pair_int_quad_map() {
chimes_ptr->build_pair_int_quad_map();
}

void chimes_compute_2b_props_fromf90(double *rij, double dr[3], char *type1, char *type2, double force[2][3], double stress[9], double *epot)
Expand Down Expand Up @@ -153,6 +163,7 @@ void chimes_compute_2b_props(double rij, double dr[3], char *atype2b[2], double
}

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) {

// convert all doubles, etc., from C to type vector for C++
vector <double> dr_3b_vec(3);
dr_3b_vec[0] = dr_3b[0];
Expand Down Expand Up @@ -190,7 +201,7 @@ void chimes_compute_3b_props(double dr_3b[3], double dist_3b[3][3], char *atype3
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) ;

Expand All @@ -215,6 +226,7 @@ void chimes_compute_3b_props(double dr_3b[3], double dist_3b[3][3], char *atype3
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)
Expand Down
2 changes: 2 additions & 0 deletions chimesFF/api/chimescalc_C.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@ int get_chimes_4b_order();
void set_chimes();
void init_chimes(int rank);
void chimes_read_params(char *param_file);
void chimes_build_pair_int_trip_map();
void chimes_build_pair_int_quad_map();
void chimes_compute_2b_props(double rij, double dr[3], char *atype2b[2], double force[2][3], double stress[9], double *epot);
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);
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);
Expand Down
8 changes: 8 additions & 0 deletions chimesFF/api/chimescalc_F.f90
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,14 @@ subroutine f_chimes_read_params(param_file) &
implicit none
character (kind=C_char), dimension(*) :: param_file
end subroutine f_chimes_read_params

subroutine f_chimes_build_pair_int_trip_map() &
& bind (C, name='chimes_build_pair_int_trip_map')
end subroutine f_chimes_build_pair_int_trip_map

subroutine f_chimes_build_pair_int_quad_map() &
& bind (C, name='chimes_build_pair_int_quad_map')
end subroutine f_chimes_build_pair_int_quad_map

function f_get_chimes_2b_order () result (order2b) &
bind (C, name='get_chimes_2b_order')
Expand Down
10 changes: 8 additions & 2 deletions chimesFF/api/chimescalc_py.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,13 @@ def read_params(param_file):
chimes_wrapper.chimes_read_params(in_param_file)
return

def build_pair_int_trip_map():
chimes_wrapper.chimes_build_pair_int_trip_map()
return

def build_pair_int_quad_map():
chimes_wrapper.chimes_build_pair_int_quad_map()
return
def chimes_compute_2b_props(rij, dr, atype2b, force, stress, epot):
"""
Compute 2-body contributions to force, stress, and energy for a
Expand Down Expand Up @@ -129,14 +135,14 @@ def chimes_compute_3b_props(dr_3b, dist_3b, atype3b, force, stress, epot):
in_dr = (ctypes.c_double * 3)(dr_3b[0], dr_3b[1], dr_3b[2])
in_dist = ((ctypes.c_double * 3) * 3)((dist_3b[0][0], dist_3b[0][1], dist_3b[0][2]),
(dist_3b[1][0], dist_3b[1][1], dist_3b[1][2]),
(dist_3b[2][0], dist_3b[2][1], dist_3b[2][2]))
(dist_3b[2][0], dist_3b[2][1], dist_3b[2][2]))
in_atype = (ctypes.c_char_p * 3)(atype3b[0].encode(), atype3b[1].encode(), atype3b[2].encode())
in_force = ((ctypes.c_double * 3) * 3) ((force[0][0], force[0][1], force[0][2]),
(force[1][0], force[1][1], force[1][2]),
(force[2][0], force[2][1], force[2][2]))
in_stress = (ctypes.c_double * 9)(stress[0], stress[1], stress[2], stress[3], stress[4], stress[5], stress[6], stress[7], stress[8])
in_epot = ctypes.c_double(epot)

chimes_wrapper.chimes_compute_3b_props(
in_dr,
in_dist,
Expand Down
4 changes: 2 additions & 2 deletions chimesFF/examples/c/README
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,9 @@ Compile with:
make all
To test the executable:

./test_wrapper-C <parameter file> <xyz file>
./C_wrapper-direct_interface <parameter file> <xyz file>

e.g. ./test_wrapper-C ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf
e.g. ./C_wrapper-direct_interface ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf


Note: xyz file must be orthorhombic, must provide lattice vectors in the comment line, e.g. lx 0.0 0.0 0.0 ly 0.0 0.0 0.0 lz, and must not
Expand Down
2 changes: 2 additions & 0 deletions chimesFF/examples/c/main.c
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ int main (int argc, char **argv)
set_chimes();
init_chimes(0);
chimes_read_params(argv[1]);
chimes_build_pair_int_trip_map();
chimes_build_pair_int_quad_map();
int order2b, order3b, order4b;
order2b = get_chimes_2b_order();
order3b = get_chimes_3b_order();
Expand Down
6 changes: 3 additions & 3 deletions chimesFF/examples/cpp/README
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,12 @@ an integer defining how many layers of ghost atoms should be used
self-interaction across the periodic boundary.

Compile with:
make test-CPP
make all

To test the executable:

./test-CPP <parameter file> <xyz file> <nlayers>
./chimescalc <parameter file> <xyz file> <nlayers>

e.g. ./test-CPP ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf 2
e.g. ./chimescalc ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf 2

Note: xyz file must be orthorhombic and must provide lattice vectors in the comment line, e.g. lx 0.0 0.0 0.0 ly 0.0 0.0 0.0 lz
2 changes: 1 addition & 1 deletion chimesFF/examples/fortran/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ chimescalc_C.o : $(CC_WRAPPER_SRC) $(CC_WRAPPER_HDR) $(CHIMESFF_SRC) $(CHIMESFF

FCC = gfortran -O3 -fPIC -std=f2003

FCC_WRAPPER_SRC=$(FCC_LOC)/../../api/chimescalc_F.F90
FCC_WRAPPER_SRC=$(FCC_LOC)/../../api/chimescalc_F.f90

chimescalc_F.o chimescalc.mod : $(FCC_WRAPPER_SRC)
$(FCC) -c $(FCC_WRAPPER_SRC) -o chimescalc_F.o
Expand Down
4 changes: 2 additions & 2 deletions chimesFF/examples/fortran/README
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ Compile with:

To test the executable:

./test_wrapper-F <parameter file> <xyz file>
./chimescalc-test_direct-F <parameter file> <xyz file>

e.g. ./test_wrapper-F ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf
e.g. ./chimescalc-test_direct-F ../../../serial_interface/tests/force_fields/published_params.liqC.2b.cubic.txt ../../../serial_interface/tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf

Note: xyz file must be orthorhombic and must provide lattice vectors in the comment line, e.g. lx 0.0 0.0 0.0 ly 0.0 0.0 0.0 lz
2 changes: 2 additions & 0 deletions chimesFF/examples/fortran/main.F90
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,8 @@ program test_F_api
call f_init_chimes(rank)
c_file = string2Cstring(param_file)
call f_chimes_read_params(c_file)
call f_chimes_build_pair_int_trip_map()
call f_chimes_build_pair_int_quad_map()
order2b=f_get_chimes_2b_order();
order3b=f_get_chimes_3b_order();
order4b=f_get_chimes_4b_order();
Expand Down
2 changes: 1 addition & 1 deletion chimesFF/examples/python/Makefile
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
CC_LOC = $(realpath .)

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

CHIMESFF_SRC=$(CC_LOC)/../../src/chimesFF.cpp
CHIMESFF_HDR=$(CC_LOC)/../../src/chimesFF.h
Expand Down
10 changes: 6 additions & 4 deletions chimesFF/examples/python/main.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def get_dist(lx,ly,lz,xrd,ycrd,zcrd,i,j):

# Initialize the ChIMES calculator

chimescalc_py.chimes_wrapper = chimescalc_py.init_chimes_wrapper("libchimescalc-direct_dl.so")
chimescalc_py.chimes_wrapper = chimescalc_py.init_chimes_wrapper(os.getcwd() + "/libchimescalc-direct_dl.so")
chimescalc_py.set_chimes()
chimescalc_py.init_chimes()

Expand All @@ -66,7 +66,9 @@ def get_dist(lx,ly,lz,xrd,ycrd,zcrd,i,j):

# Read the parameters

wrapper_py.read_params(param_file)
chimescalc_py.read_params(param_file)
chimescalc_py.build_pair_int_trip_map()
chimescalc_py.build_pair_int_quad_map()

# Read the coordinates, set up the force, stress, and energy vars

Expand Down Expand Up @@ -129,7 +131,7 @@ def get_dist(lx,ly,lz,xrd,ycrd,zcrd,i,j):

forces[i][0] = tmp_force[0][0]; forces[i][1] = tmp_force[0][1]; forces[i][2] = tmp_force[0][2]
forces[j][0] = tmp_force[1][0]; forces[j][1] = tmp_force[1][1]; forces[j][2] = tmp_force[1][2]
if (order_3b > 0) or (order_4b> 0):

for k in range(j+1,natoms):
Expand All @@ -153,7 +155,7 @@ def get_dist(lx,ly,lz,xrd,ycrd,zcrd,i,j):
forces[i][0] = tmp_force[0][0]; forces[i][1] = tmp_force[0][1]; forces[i][2] = tmp_force[0][2]
forces[j][0] = tmp_force[1][0]; forces[j][1] = tmp_force[1][1]; forces[j][2] = tmp_force[1][2]
forces[k][0] = tmp_force[2][0]; forces[k][1] = tmp_force[2][1]; forces[k][2] = tmp_force[2][2]

if dist_ik >= maxcut_4b:
continue
if dist_jk >= maxcut_4b:
Expand Down
13 changes: 7 additions & 6 deletions chimesFF/src/chimesFF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1562,23 +1562,23 @@ void chimesFF::compute_3B(const vector<double> & dx, const vector<double> & dr,
}
#endif


int type_idx = typ_idxs[0]*natmtyps*natmtyps + typ_idxs[1]*natmtyps + typ_idxs[2] ;
int tripidx = atom_int_trip_map[type_idx];

if(tripidx < 0) // Skipping an excluded interaction
return;

// Check whether cutoffs are within allowed ranges
vector<int> & mapped_pair_idx = pair_int_trip_map[type_idx] ;



if (dx[0] >= chimes_3b_cutoff[ tripidx ][1][mapped_pair_idx[0]]) // ij
return;
if (dx[1] >= chimes_3b_cutoff[ tripidx ][1][mapped_pair_idx[1]]) // ik
return;
if (dx[2] >= chimes_3b_cutoff[ tripidx ][1][mapped_pair_idx[2]]) // jk
return;

// At this point, all distances are within allowed ranges. We can now proceed to the force/stress/energy calculation

#ifdef USE_DISTANCE_TENSOR
Expand All @@ -1587,6 +1587,7 @@ void chimesFF::compute_3B(const vector<double> & dx, const vector<double> & dr,
init_distance_tensor(dr2, dr, npairs) ;
#endif


// Set up the polynomials

set_cheby_polys(Tn_ij, Tnd_ij, dx[0], atom_int_pair_map[ typ_idxs[0]*natmtyps + typ_idxs[1] ], chimes_3b_cutoff[tripidx][0][mapped_pair_idx[0]], chimes_3b_cutoff[tripidx][1][mapped_pair_idx[0]], 1);
Expand All @@ -1611,7 +1612,7 @@ void chimesFF::compute_3B(const vector<double> & dx, const vector<double> & dr,
double coeff;
int powers[npairs] ;
double force_scalar[npairs] ;

for(int coeffs=0; coeffs<ncoeffs_3b[tripidx]; coeffs++)
{
coeff = chimes_3b_params[tripidx][coeffs];
Expand Down Expand Up @@ -1710,7 +1711,7 @@ void chimesFF::compute_3B(const vector<double> & dx, const vector<double> & dr,
stress[5] -= force_scalar[2] * dr[2*CHDIM+2] * dr[2*CHDIM+2]; // zz tensor component
#endif
}

force_scalar_in[0] = force_scalar[0];
force_scalar_in[1] = force_scalar[1];
force_scalar_in[2] = force_scalar[2];
Expand Down
34 changes: 27 additions & 7 deletions chimesFF/src/chimesFF.h
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,6 @@ class chimesFF
void compute_4B(const vector<double> & dx, const vector<double> & dr, const vector<int> & typ_idxs, vector<double> & force, vector<double> & stress, double & energy, chimes4BTmp &tmp);
void compute_4B(const vector<double> & dx, const vector<double> & dr, const vector<int> & typ_idxs, vector<double> & force, vector<double> & stress, double & energy, chimes4BTmp &tmp, vector<double> & force_scalar_in);


void get_cutoff_2B(vector<vector<double> > & cutoff_2b); // Populates the 2b cutoffs

double max_cutoff_2B(bool silent = false); // Returns the largest 2B cutoff
Expand All @@ -175,6 +174,11 @@ class chimesFF
void build_pair_int_trip_map() ;
void build_pair_int_quad_map() ;

// Functions to aid using ChIMES Calculator for fitting

inline int get_badness();
inline void reset_badness();

private:

string xform_style; // Morse, direct, inverse, etc...
Expand All @@ -184,6 +188,7 @@ class chimesFF
vector<double> morse_var; // [npairs]; morse_lambda
vector<double> penalty_params; // [2]; Second dimension: [0] = A_pen, [1] = d_pen
vector<double> energy_offsets; // [natmtyps]; Single atom ChIMES energies
int badness; // Keeps track of whether any interactions for atoms owned by proc rank are below rcutin, in the penalty region, or in the r>rcutin+dp region. 0 = good, 1 = in penalty region, 2 = below rcutin

// Names (chemical symbols for constituent atoms) .. handled differently for 2-body versus >2-body interactions

Expand Down Expand Up @@ -339,28 +344,43 @@ inline void chimesFF::get_penalty(const double dx, const int & pair_idx, double

E_penalty = 0.0;
force_scalar = 1.0;
if (dx - penalty_params[0] < chimes_2b_cutoff[pair_idx][0])


if (dx - penalty_params[0] < chimes_2b_cutoff[pair_idx][0]) // Then we're within the penalty-enforced region of distance space
{
r_penalty = chimes_2b_cutoff[pair_idx][0] + penalty_params[0] - dx;

if(dx < chimes_2b_cutoff[pair_idx][0])
badness = 2;
else if (1 > badness) // Only update badness if candiate badness is worse than its current value
badness = 1;
}
if ( r_penalty > 0.0 )
{
E_penalty = r_penalty * r_penalty * r_penalty * penalty_params[1];

force_scalar = -3.0 * r_penalty * r_penalty * penalty_params[1];

if (rank == 0)
{
//if (rank == 0) // Commenting out - we need all ranks to report if the penalty function has been sampled
//{
cout << "chimesFF: " << "Adding penalty in 2B Cheby calc, r < rmin+penalty_dist " << fixed
<< dx << " "
<< chimes_2b_cutoff[pair_idx][0] + penalty_params[0]
<< " pair type: " << pair_idx << endl;
cout << "chimesFF: " << "\t...Penalty potential = "<< E_penalty << endl;
}
//}
}
}

inline int chimesFF::get_badness()
{
return badness;
}

inline void chimesFF::reset_badness()
{
badness = 0;
}

inline void chimesFF::build_atom_and_pair_mappers(const int natoms, const int npairs, const vector<int> & typ_idxs,
const vector<string> & clu_params_pair_typs, vector<int> & mapped_pair_idx)
// Interface to array-based version.
Expand Down
Loading

0 comments on commit 6cb12c0

Please sign in to comment.