Printing Pairwise SP EFP Energies

.gitignore

@@ -19,3 +19,5 @@
 /efpmd/src/efpmd
 /fraglib/params
 tags
+bin/*
+include/*

efpmd/src/Makefile

@@ -6,7 +6,7 @@ LIBS= -lefp -lopt -lff $(MYLIBS) -lm
 
 PROG= efpmd
 ALL_O= cfg.o common.o efield.o energy.o grad.o gtest.o hess.o main.o \
-       md.o msg.o opt.o parse.o rand.o sp.o
+       md.o msg.o opt.o parse.o rand.o sp.o psp.o mc_direct.o
 
 $(PROG): $(ALL_O)
 	$(CC) -o $@ $(CFLAGS) $(LDFLAGS) $(ALL_O) $(LIBS)

efpmd/src/common.c

@@ -119,13 +119,13 @@ void print_geometry(struct efp *efp)
 		}
 	}
 
+
 	size_t n_charges;
 	check_fail(efp_get_point_charge_count(efp, &n_charges));
-
+ 
 	if (n_charges > 0) {
 		double xyz[3 * n_charges];
 		check_fail(efp_get_point_charge_coordinates(efp, xyz));
-
 		for (size_t i = 0; i < n_charges; i++) {
 			char label[32];
 			double x = xyz[3 * i + 0] * BOHR_RADIUS;
@@ -136,7 +136,7 @@ void print_geometry(struct efp *efp)
 			msg("%-16s %12.6lf %12.6lf %12.6lf\n", label, x, y, z);
 		}
 	}
-
+	
 	msg("\n\n");
 }
 
@@ -280,3 +280,7 @@ vec_t box_from_str(const char *str)
 	vec_scale(&box, 1.0 / BOHR_RADIUS);
 	return box;
 }
+
+void print_six_t(const double *six_t){
+         msg("   %16.8E %16.8E %16.8E %16.8E %16.8E %16.8E", six_t[1], six_t[3], six_t[4], six_t[2], six_t[5], six_t[0]);
+}

efpmd/src/common.h

@@ -59,7 +59,9 @@ enum run_type {
 	RUN_TYPE_OPT,
 	RUN_TYPE_MD,
 	RUN_TYPE_EFIELD,
-	RUN_TYPE_GTEST
+	RUN_TYPE_GTEST,
+	RUN_TYPE_PSP,
+	RUN_TYPE_MC
 };
 
 enum ensemble_type {
@@ -96,6 +98,7 @@ struct state {
 	struct sys *sys;
 	double energy;
 	double *grad;
+	double *energy_components; 
 };
 
 void NORETURN die(const char *, ...);
@@ -113,6 +116,7 @@ void print_fragment(const char *, const double *, const double *);
 void print_charge(double, double, double, double);
 void print_vector(size_t, const double *);
 void print_matrix(size_t, size_t, const double *);
+void print_six_t(const double *); 
 
 void check_fail(enum efp_result);
 void compute_energy(struct state *, bool);

efpmd/src/energy.c

@@ -36,6 +36,7 @@ void compute_energy(struct state *state, bool do_grad)
 	int itotal;
 
 	/* EFP part */
+
 	check_fail(efp_compute(state->efp, do_grad));
 	check_fail(efp_get_energy(state->efp, &efp_energy));
 	check_fail(efp_get_frag_count(state->efp, &nfrag));
@@ -46,6 +47,10 @@ void compute_energy(struct state *state, bool do_grad)
 		    state->grad + 6 * nfrag));
 	}
 
+        if (cfg_get_bool(state->cfg, "enable_pairwise")){
+                check_fail(efp_get_energy_components(state->efp, state->energy_components));
+        }
+
 	state->energy = efp_energy.total;
 
 	/* constraints */

efpmd/src/main.c

@@ -37,6 +37,8 @@ void sim_opt(struct state *);
 void sim_md(struct state *);
 void sim_efield(struct state *);
 void sim_gtest(struct state *);
+void sim_psp(struct state *);
+void sim_mc(struct state *);
 
 #define USAGE_STRING \
 	"usage: efpmd [-d | -v | -h | input]\n" \
@@ -55,14 +57,18 @@ static struct cfg *make_cfg(void)
 		"opt\n"
 		"md\n"
 		"efield\n"
-		"gtest\n",
+		"gtest\n"
+		"psp\n"
+		"mc\n",
 		(int []) { RUN_TYPE_SP,
 			   RUN_TYPE_GRAD,
 			   RUN_TYPE_HESS,
 			   RUN_TYPE_OPT,
 			   RUN_TYPE_MD,
 			   RUN_TYPE_EFIELD,
-			   RUN_TYPE_GTEST });
+			   RUN_TYPE_GTEST,
+			   RUN_TYPE_PSP,
+			   RUN_TYPE_MC });
 
 	cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
 		"xyzabc\n"
@@ -138,7 +144,11 @@ static struct cfg *make_cfg(void)
 	cfg_add_double(cfg, "pressure", 1.0);
 	cfg_add_double(cfg, "thermostat_tau", 1.0e3);
 	cfg_add_double(cfg, "barostat_tau", 1.0e4);
-
+	cfg_add_int(cfg, "ligand", 0); 
+	cfg_add_bool(cfg, "enable_pairwise", false); 
+	cfg_add_double(cfg, "dispmag_threshold", 0.5); 
+	cfg_add_double(cfg, "dispmag_modifier", 0.9); 
+	cfg_add_int(cfg, "dispmag_modify_steps", 100); 
 	return cfg;
 }
 
@@ -159,6 +169,11 @@ static sim_fn_t get_sim_fn(enum run_type run_type)
 		return sim_efield;
 	case RUN_TYPE_GTEST:
 		return sim_gtest;
+    case RUN_TYPE_PSP:
+        return sim_psp;
+    case RUN_TYPE_MC:
+    	fprintf(stdout, "\n----monte carlo----\n\n");
+        return sim_mc;
 	}
 	assert(0);
 }
@@ -247,7 +262,9 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 		.pol_driver = cfg_get_enum(cfg, "pol_driver"),
 		.enable_pbc = cfg_get_bool(cfg, "enable_pbc"),
 		.enable_cutoff = cfg_get_bool(cfg, "enable_cutoff"),
-		.swf_cutoff = cfg_get_double(cfg, "swf_cutoff")
+		.swf_cutoff = cfg_get_double(cfg, "swf_cutoff"),
+		.enable_pairwise = cfg_get_bool(cfg, "enable_pairwise"), 
+		.ligand = cfg_get_int(cfg, "ligand")
 	};
 
 	enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord");
@@ -309,6 +326,7 @@ static void state_init(struct state *state, const struct cfg *cfg, const struct
 	state->energy = 0;
 	state->grad = xcalloc(sys->n_frags * 6 + sys->n_charges * 3, sizeof(double));
 	state->ff = NULL;
+	state->energy_components = xcalloc(sys->n_frags * 6, sizeof(double));
 
 	if (cfg_get_bool(cfg, "enable_ff")) {
 		if ((state->ff = ff_create()) == NULL)