Merge pull request #57 from rk-lindsey/develop

Update main from develop

CMakeLists.txt

@@ -92,6 +92,11 @@ target_link_libraries     (chimescalc-test_serial-C ChimesCalc)
 target_compile_features   (chimescalc-test_serial-C PRIVATE cxx_std_11)
 target_compile_definitions(chimescalc-test_serial-C PRIVATE "DEBUG=${DEBUG}")
 
+add_executable(chimescalc-test_serial-C_instance serial_interface/examples/c_instance/main.c)
+target_link_libraries     (chimescalc-test_serial-C_instance ChimesCalc)
+target_compile_features   (chimescalc-test_serial-C_instance PRIVATE cxx_std_11)
+target_compile_definitions(chimescalc-test_serial-C_instance PRIVATE "DEBUG=${DEBUG}")
+
 
 ####################################################################################################
 # Dynamically loadable library (e.g for Python)
@@ -163,6 +168,9 @@ if(WITH_FORTRAN_API OR WITH_FORTRAN08_API)
     add_executable(chimescalc-test_serial-F serial_interface/examples/fortran/main.F90)
     target_link_libraries(chimescalc-test_serial-F ChimesCalc_Fortran)
     target_compile_definitions(chimescalc-test_serial-F PRIVATE "DEBUG=${DEBUG}")
+    add_executable(chimescalc-test_serial-F_instance serial_interface/examples/fortran_instance/main.F90)
+    target_link_libraries(chimescalc-test_serial-F_instance ChimesCalc_Fortran)
+    target_compile_definitions(chimescalc-test_serial-F_instance PRIVATE "DEBUG=${DEBUG}")
 
     # Optional, as ancient Fortran compilers may have difficulties with it
     if(WITH_FORTRAN08_API)
@@ -210,7 +218,9 @@ enable_testing()
 set(apis
     "cpp;${CMAKE_CURRENT_BINARY_DIR}/chimescalc"
     "c;${CMAKE_CURRENT_BINARY_DIR}/chimescalc-test_serial-C"
+    "c_instance;${CMAKE_CURRENT_BINARY_DIR}/chimescalc-test_serial-C_instance"
     "fortran;${CMAKE_CURRENT_BINARY_DIR}/chimescalc-test_serial-F"
+    "fortran_instance;${CMAKE_CURRENT_BINARY_DIR}/chimescalc-test_serial-F_instance"
     "fortran08;${CMAKE_CURRENT_BINARY_DIR}/chimescalc-test_serial-F08")
 
 set(_testdir "${CMAKE_CURRENT_SOURCE_DIR}/serial_interface/tests")

CODEOWNERS

@@ -0,0 +1,2 @@
+* @RKLindsey
+* @rk-lindsey

chimesFF/api/chimescalc_C.cpp

@@ -71,6 +71,8 @@ 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() {

chimesFF/api/chimescalc_py.py

@@ -135,7 +135,7 @@ 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]),

chimesFF/src/chimesFF.h

@@ -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
@@ -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...
@@ -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
 
@@ -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.

etc/lmp/src/pair_chimes.cpp

@@ -34,10 +34,12 @@
 #include "error.h"
 #include "pair_chimes.h"
 #include "group.h"
-#include "update.h"
-
+#include "update.h" // Needed for mb neighlist updates and info printing for fitting
+#include "output.h" // Needed for infor printing for fitting -- dump 1 must be the "main" dump file used for fitting
+#include  "utils.h" // Needed for infor printing for fitting
 #include <vector>
 #include <iostream>
+#include <sstream>
 
 using namespace LAMMPS_NS;
 
@@ -54,9 +56,9 @@ init_style 	(done)		initialization specific to this pair style
 
 write_restart			write i,j pair coeffs to restart file
 read_restart			read i,j pair coeffs from restart file
-write_restart_settings		write global settings to restart file
-read_restart_settings		read global settings from restart file
-single				force and energy fo a single pairwise interaction between two atoms
+write_restart_settings	write global settings to restart file
+read_restart_settings	read global settings from restart file
+single				    force and energy fo a single pairwise interaction between two atoms
 */
 
 
@@ -67,7 +69,8 @@ PairCHIMES::PairCHIMES(LAMMPS *lmp) : Pair(lmp)
 	int me = comm->me;
 	MPI_Comm_rank(world,&me);
 
-	chimes_calculator.init(me);
+	chimes_calculator.init(me);  
+    for_fitting = false;
 
 	// 2, 3, and 4-body vars for chimesFF access
 
@@ -110,12 +113,26 @@ PairCHIMES::~PairCHIMES()
 	    memory->destroy(setflag);
 	    memory->destroy(cutsq);
 	}
+
+    if (badness_stream.is_open())
+        badness_stream.close();
 }	
 
 void PairCHIMES::settings(int narg, char **arg)
 {
-	if (narg != 0) 
-		error -> all(FLERR,"Illegal pair_style command. Expects no arguments beyond pair_style name.");
+	if (narg > 1) 
+		error -> all(FLERR,"Illegal pair_style command. Expects no more than one arguments (string: fitting)");
+
+    if (narg == 1) 
+    {  
+        if (utils::strmatch(arg[0],"fitting"))
+        {
+            for_fitting   = true;
+            stringstream ss;
+            ss << chimes_calculator.rank;
+            badness_stream.open("rank-" + ss.str() + ".badness.log");    
+        }
+    }
 
 	return;	
 }
@@ -132,6 +149,7 @@ void PairCHIMES::coeff(int narg, char **arg)
 	chimes_calculator.read_parameters(chimesFF_paramfile);
 
 	set_chimes_type();
+
     //chimes_calculator.set_atomtypes(chimes_type);
     chimes_calculator.build_pair_int_trip_map() ; 
     chimes_calculator.build_pair_int_quad_map() ;
@@ -238,11 +256,11 @@ void PairCHIMES::build_mb_neighlists()
 	neighborlist_3mers.clear();
 	neighborlist_4mers.clear();
 
-	int i,j,k,l,inum,jnum,knum,lnum, ii, jj, kk, ll;		// Local iterator vars
-	int *ilist,*jlist,*klist,*llist, *numneigh,**firstneigh;	// Local neighborlist vars
-	tagint 	*tag   = atom -> tag;					// Access to global atom indices
-	int     itag, jtag, ktag, ltag;					// holds tags	
-	double 	**x    = atom -> x;					// Access to system coordinates
+	int i,j,k,l,inum,jnum,knum,lnum, ii, jj, kk, ll;		 // Local iterator vars
+	int *ilist,*jlist,*klist,*llist, *numneigh,**firstneigh; // Local neighborlist vars
+	tagint 	*tag   = atom -> tag;					         // Access to global atom indices
+	int     itag, jtag, ktag, ltag;					         // holds tags	
+	double 	**x    = atom -> x;					             // Access to system coordinates
 
 	double maxcut_3b_padded = maxcut_3b + neighbor-> skin;
 	double maxcut_4b_padded = maxcut_4b + neighbor-> skin;
@@ -253,12 +271,12 @@ void PairCHIMES::build_mb_neighlists()
 	// Access to neighbor list vars
 	////////////////////////////////////////
 
-	inum       = list -> inum; 		// length of the list
-	ilist      = list -> ilist; 		// list of i atoms for which neighbor list exists
-	numneigh   = list -> numneigh;		// length of each of the ilist neighbor lists
-	firstneigh = list -> firstneigh;	// point to the list of neighbors of i	
+	inum       = list -> inum; 		 // length of the list
+	ilist      = list -> ilist; 	 // list of i atoms for which neighbor list exists
+	numneigh   = list -> numneigh;	 // length of each of the ilist neighbor lists
+	firstneigh = list -> firstneigh; // point to the list of neighbors of i	
 
-	for (ii = 0; ii < inum; ii++)	// Loop over real atoms (ai)	
+	for (ii = 0; ii < inum; ii++) // Loop over real atoms (ai)	
 	{
 		i     = ilist[ii];		
 		itag  = tag[i];			
@@ -385,26 +403,26 @@ void PairCHIMES::compute(int eflag, int vflag)
 
 	// Temp vars to hold chimes output for passing to ev_tally function
 
-    std::vector<double>         fscalar(6);
-    std::vector<double>         tmp_dist(1);
-    std::vector<double>         tmp_dr(6);
-    int atmidxlst[6][2];
+    std::vector<double>  fscalar(6);
+    std::vector<double>  tmp_dist(1);
+    std::vector<double>  tmp_dr(6);
+    int                  atmidxlst[6][2];
 
 	// General LAMMPS compute vars
 
-	int 	i,j,k,l,inum,jnum, ii, jj;		// Local iterator vars
+	int 	i,j,k,l,inum,jnum, ii, jj;	// Local iterator vars
 	int 	*ilist,*jlist, *numneigh,**firstneigh;	// Local neighborlist vars
 	int     idx;
 
-	double 	**x    = atom -> x;		// Access to system coordinates
-	double 	**f    = atom -> f;		// Access to system forces
+	double 	**x    = atom -> x;		    // Access to system coordinates
+	double 	**f    = atom -> f;		    // Access to system forces
 
 	int 	*type  = atom -> type;		// Acces to system atom types (countng starts from 1, chimesFF class expects counting from 0!)
 	tagint 	*tag   = atom -> tag;		// Access to global atom indices (sort of like "parent" indices)
 	int     itag, jtag, ktag, ltag;		// holds tags
 	int 	nlocal = atom -> nlocal;	// Number of real atoms owned by current process .. used used to assure force assignments aren't duplicated
 	int 	newton_pair = force -> newton_pair;	// Should f_j be automatically set to -f_i (true) or manually calculated (false)
-	double  energy;				// pair energy 
+	double  energy;				        // pair energy 
 
 	int me = comm->me;
 	MPI_Comm_rank(world,&me);	
@@ -426,10 +444,10 @@ void PairCHIMES::compute(int eflag, int vflag)
 	// Access to (2-body) neighbor list vars
 	////////////////////////////////////////
 
-	inum       = list -> inum; 		// length of the list
-	ilist      = list -> ilist; 		// list of i atoms for which neighbor list exists
-	numneigh   = list -> numneigh;		// length of each of the ilist neighbor lists
-	firstneigh = list -> firstneigh;	// point to the list of neighbors of i
+	inum       = list -> inum; 		 // length of the list
+	ilist      = list -> ilist; 	 // list of i atoms for which neighbor list exists
+	numneigh   = list -> numneigh;	 // length of each of the ilist neighbor lists
+	firstneigh = list -> firstneigh; // point to the list of neighbors of i
 
 
     chimes2BTmp chimes_2btmp(chimes_calculator.poly_orders[0]) ;
@@ -451,6 +469,10 @@ void PairCHIMES::compute(int eflag, int vflag)
 			std::cout << "	...update complete" << std::endl;
 		}
 	}
+
+    // Prepare the badness variable
+
+    chimes_calculator.reset_badness();
 
 	////////////////////////////////////////
 	// Compute 1- and 2-body interactions
@@ -475,7 +497,7 @@ void PairCHIMES::compute(int eflag, int vflag)
 
 		// Now move on to two-body force, stress, and energy
 
-		for (jj = 0; jj < jnum; jj++)			// Loop over neighbors of i
+		for (jj = 0; jj < jnum; jj++) // Loop over neighbors of i
 		{
 			j     = jlist[jj];			// Index of the jj atom
 			jtag  = tag[j];				// Get j's global atom index (sort of like its "parent")
@@ -489,7 +511,7 @@ void PairCHIMES::compute(int eflag, int vflag)
 
 			dist = get_dist(i,j,&dr[0]);
 
-			typ_idxs_2b[0] = chimes_type[type[i]-1];		// Type (index) of the current atom... subtract 1 to account for chimesFF vs LAMMPS numbering convention
+			typ_idxs_2b[0] = chimes_type[type[i]-1]; // Type (index) of the current atom... subtract 1 to account for chimesFF vs LAMMPS numbering convention
 			typ_idxs_2b[1] = chimes_type[type[j]-1];
 
 			// Using std::fill for maximum efficiency.
@@ -500,7 +522,7 @@ void PairCHIMES::compute(int eflag, int vflag)
 
 			energy = 0.0;	
 
-			chimes_calculator.compute_2B( dist, dr, typ_idxs_2b, force_2b, stensor, energy, chimes_2btmp);			
+			chimes_calculator.compute_2B( dist, dr, typ_idxs_2b, force_2b, stensor, energy, chimes_2btmp);	// Auto-updates badness		
 
 			for (idx=0; idx<3; idx++)
 			{
@@ -522,6 +544,11 @@ void PairCHIMES::compute(int eflag, int vflag)
 				ev_tally_mb(2, atmidxlst, energy, fscalar, tmp_dist, dr);         
 		}
 	}
+
+    // Document badness for configuration: current timestep, current rank, worst badness seen by rank
+    if (for_fitting)
+        if(update->ntimestep % output->every_dump[0] == 0)
+            badness_stream << update->ntimestep << " " <<  chimes_calculator.get_badness() << endl;
 
 	if (chimes_calculator.poly_orders[1] > 0)
 	{
@@ -644,11 +671,11 @@ void PairCHIMES::set_chimes_type()
 {
 	int nmatches = 0;
 
-	for (int i=1; i<= atom->ntypes; i++) 						// Lammps indexing starts at 1
+	for (int i=1; i<= atom->ntypes; i++) // Lammps indexing starts at 1
 	{
-		for (int j=0; j<chimes_calculator.natmtyps; j++) 			// ChIMES indexing starts at 0
+		for (int j=0; j<chimes_calculator.natmtyps; j++) // ChIMES indexing starts at 0
 		{
-			if (abs(atom->mass[i] - chimes_calculator.masses[j]) < 1e-3)	// Masses should match to at least 3 decimal places
+			if (abs(atom->mass[i] - chimes_calculator.masses[j]) < 1e-3) // Masses should match to at least 3 decimal places
 			{
 				chimes_type.push_back(j);
 				nmatches++;

etc/lmp/src/pair_chimes.h

@@ -42,9 +42,9 @@ init_style 	(done)		initialization specific to this pair style
 
 write_restart			write i,j pair coeffs to restart file
 read_restart			read i,j pair coeffs from restart file
-write_restart_settings		write global settings to restart file
-read_restart_settings		read global settings from restart file
-single				force and energy fo a single pairwise interaction between two atoms
+write_restart_settings	write global settings to restart file
+read_restart_settings	read global settings from restart file
+single				    force and energy fo a single pairwise interaction between two atoms
 */
 
 
@@ -56,21 +56,25 @@ namespace LAMMPS_NS
 
 			// Variable definitions
 
-			chimesFF chimes_calculator; // chimesFF instance
+			chimesFF chimes_calculator;   // chimesFF instance
 
-			char * chimesFF_paramfile;	// ChIMES parameter file
+			char * chimesFF_paramfile;	  // ChIMES parameter file
 
-			std::vector<int> chimes_type;	// For i = LMP atom type indx, chimes_type[i-1] gives the ChIMES parameter file type idx
+			std::vector<int> chimes_type; // For i = LMP atom type indx, chimes_type[i-1] gives the ChIMES parameter file type idx
 
 			double maxcut_3b;
 			double maxcut_4b;
 
-			int n_3mers;				// number of neighborlist_Xmers entries
+			int n_3mers;				   // number of neighborlist_Xmers entries
 			int n_4mers;
 
 			std::vector<std::vector<int> > neighborlist_3mers;	// custom neighbor list; neighborlist_Xmers[cluster idx][atom in cluster idx]
 			std::vector<std::vector<int> > neighborlist_4mers;
-
+
+            // Prepare files necessary for ChIMES fitting 
+
+            bool     for_fitting;
+            ofstream badness_stream;			
 
 			// 2-body vars for chimesFF access