conflict resolution

Author: slipchenko

.gitignore

@@ -13,9 +13,11 @@
 .deps/
 .dirstamp
 .libs/
+.idea
 /config.inc
 /doxygen_html
 /efpmd/libff/libtinker.a
 /efpmd/src/efpmd
 /fraglib/params
+/bin/efpmd
 tags

.idea/libefp.iml

@@ -215,6 +215,20 @@ Default value: `"."` (current directory)
 The `<path>` parameter should not contain spaces or should be in double quotes
 otherwise.
 
+### Pairwise energy analysis
+
+##### Enable/Disable pairwise analysis
+
+`enable_pairwise [true|false]`
+
+Default value: `false`
+
+##### Specify ligand 
+
+`ligand [integer]`
+
+Default value: `0` (the first fragment in the system)
+
 ### Periodic Boundary Conditions (PBC)
 
 ##### Enable/Disable PBC

efpmd/README.md

@@ -168,6 +168,7 @@ void print_energy(struct state *state)
 	msg("\n\n");
 }
 
+
 void print_gradient(struct state *state)
 {
 	size_t n_frags;
@@ -270,6 +271,98 @@ void print_matrix(size_t rows, size_t cols, const double *mat)
 	}
 }
 
+void print_pair_energy(struct state *state){
+
+	if (cfg_get_bool(state->cfg, "enable_pairwise")) {
+		msg(" ------ PAIRWISE ENERGY ANALYSIS FOLLOWS ------------------\n\n");
+		size_t n_frags;
+        check_fail(efp_get_frag_count(state->efp, &n_frags));
+        double coord[6 * n_frags];
+        check_fail(efp_get_coordinates(state->efp, coord));
+
+        struct efp_energy *energies;
+        energies = xmalloc(n_frags * sizeof(struct efp_energy));
+        check_fail(efp_get_pairwise_energy(state->efp, energies));
+
+        char ligand[32];
+        size_t lig_atoms;
+
+        size_t ligand_index = cfg_get_int(state->cfg, "ligand");
+        check_fail(efp_get_frag_name(state->efp, ligand_index, sizeof(ligand),ligand));
+        check_fail(efp_get_frag_atom_count(state->efp, ligand_index, &lig_atoms));
+        struct efp_atom latoms[lig_atoms];
+        check_fail(efp_get_frag_atoms(state->efp, ligand_index, lig_atoms, latoms));
+
+        char frag_name[32];
+        size_t frag_atoms;
+        double lattice_energy[6];
+        for (size_t j=0; j<6; j++){
+        	lattice_energy[j]=0.0;
+        }
+        for (size_t i=0; i <n_frags; i++){
+                check_fail(efp_get_frag_name(state->efp, i, sizeof(frag_name),frag_name));
+                check_fail(efp_get_frag_atom_count(state->efp, i, &frag_atoms));
+
+                struct efp_atom atoms[frag_atoms];
+                check_fail(efp_get_frag_atoms(state->efp, i, frag_atoms, atoms));
+
+                msg("   PAIRWISE ENERGY BETWEEN FRAGMENT %zu (%s) AND LIGAND %zu (%s) \n", i, frag_name, ligand_index, ligand);
+                msg("fragment %s\n", frag_name);
+                for (size_t a = 0; a < frag_atoms; a++) {
+                        double x = atoms[a].x * BOHR_RADIUS;
+                        double y = atoms[a].y * BOHR_RADIUS;
+                        double z = atoms[a].z * BOHR_RADIUS;
+                        msg("   %-16s %12.6lf %12.6lf %12.6lf\n", atoms[a].label, x, y, z);
+                }
+                msg("\n");
+
+                msg("fragment %s\n", ligand);
+                for (size_t a = 0; a < lig_atoms; a++) {
+                        double x = latoms[a].x * BOHR_RADIUS;
+                        double y = latoms[a].y * BOHR_RADIUS;
+                        double z = latoms[a].z * BOHR_RADIUS;
+                        msg("   %-16s %12.6lf %12.6lf %12.6lf\n", latoms[a].label, x, y, z);
+                }
+                msg("\n");
+
+                msg("    PAIRWISE ENERGY COMPONENTS (ATOMIC UNITS)\n");
+                msg("%40s %16.10lf\n", "PAIRWISE ELECTROSTATIC ENERGY", energies[i].electrostatic);
+                msg("%40s %16.10lf\n", "PAIRWISE POLARIZATION ENERGY", energies[i].polarization);
+                msg("%40s %16.10lf\n", "PAIRWISE DISPERSION ENERGY", energies[i].dispersion);
+                msg("%40s %16.10lf\n", "PAIRWISE EXCHANGE REPULSION ENERGY",
+                        energies[i].exchange_repulsion);
+                msg("%40s %16.10lf\n", "PAIRWISE CHARGE PENETRATION ENERGY",
+                		energies[i].charge_penetration);
+
+                energies[i].total = energies[i].electrostatic + energies[i].polarization + energies[i].dispersion
+                		+ energies[i].exchange_repulsion + energies[i].charge_penetration;
+                msg("%40s %16.10lf\n", "PAIRWISE TOTAL ENERGY", energies[i].total);
+                msg("\n ---------------------------------------------------------\n");
+
+                lattice_energy[0] = lattice_energy[0] + energies[i].electrostatic;
+                lattice_energy[1] = lattice_energy[1] + energies[i].polarization;
+                lattice_energy[2] = lattice_energy[2] + energies[i].dispersion;
+                lattice_energy[3] = lattice_energy[3] + energies[i].exchange_repulsion;
+                lattice_energy[4] = lattice_energy[4] + energies[i].charge_penetration;
+                lattice_energy[5] = lattice_energy[5] + energies[i].total;
+        }
+        free(energies);
+
+        msg("    LATTICE ENERGY COMPONENTS (ATOMIC UNITS)\n");
+        msg("%40s %16.10lf\n", "LATTICE ELECTROSTATIC ENERGY", lattice_energy[0]);
+        msg("%40s %16.10lf\n", "LATTICE POLARIZATION ENERGY", lattice_energy[1]);
+        msg("%40s %16.10lf\n", "LATTICE DISPERSION ENERGY", lattice_energy[2]);
+        msg("%40s %16.10lf\n", "LATTICE EXCHANGE REPULSION ENERGY", lattice_energy[3]);
+        msg("%40s %16.10lf\n", "LATTICE CHARGE PENETRATION ENERGY", lattice_energy[4]);
+        msg("%40s %16.10lf\n", "LATTICE TOTAL ENERGY", lattice_energy[5]);
+        msg("\n\n");
+
+        msg(" ------ PAIRWISE ENERGY ANALYSIS COMPLETED --------------\n\n");
+	}
+}
+
+
+
 six_t box_from_str(const char *str)
 {
 	six_t box;
@@ -283,3 +376,4 @@ six_t box_from_str(const char *str)
 	// vec_scale(&box, 1.0 / BOHR_RADIUS);
 	return box;
 }
+

efpmd/src/common.c

@@ -59,8 +59,9 @@ enum run_type {
 	RUN_TYPE_OPT,
 	RUN_TYPE_MD,
 	RUN_TYPE_EFIELD,
-	RUN_TYPE_GTEST
+	RUN_TYPE_GTEST,
 };
+//RUN_TYPE_PAIRWISE
 
 enum ensemble_type {
 	ENSEMBLE_TYPE_NVE,
@@ -113,6 +114,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_pair_energy(struct state *);
 
 void check_fail(enum efp_result);
 void compute_energy(struct state *, bool);
@@ -121,4 +123,5 @@ six_t box_from_str(const char *);
 int efp_strcasecmp(const char *, const char *);
 int efp_strncasecmp(const char *, const char *, size_t);
 
+
 #endif /* EFPMD_COMMON_H */

efpmd/src/common.h

@@ -37,6 +37,7 @@ void sim_opt(struct state *);
 void sim_md(struct state *);
 void sim_efield(struct state *);
 void sim_gtest(struct state *);
+// void sim_pairwise(struct state *);
 
 #define USAGE_STRING \
 	"usage: efpmd [-d | -v | -h | input]\n" \
@@ -62,7 +63,9 @@ static struct cfg *make_cfg(void)
 			   RUN_TYPE_OPT,
 			   RUN_TYPE_MD,
 			   RUN_TYPE_EFIELD,
-			   RUN_TYPE_GTEST });
+			   RUN_TYPE_GTEST});
+	/* 		"pairwise\n",
+	 RUN_TYPE_PAIRWISE */
 
 	cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
 		"xyzabc\n"
@@ -139,6 +142,10 @@ static struct cfg *make_cfg(void)
 	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); 
+
+
 	return cfg;
 }
 
@@ -160,6 +167,8 @@ static sim_fn_t get_sim_fn(enum run_type run_type)
 	case RUN_TYPE_GTEST:
 		return sim_gtest;
 	}
+	// 	case RUN_TYPE_PAIRWISE:
+	//  return sim_ pairwise;
 	assert(0);
 }
 
@@ -247,7 +256,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");
@@ -310,6 +321,7 @@ static void state_init(struct state *state, const struct cfg *cfg, const struct
 	state->grad = xcalloc(sys->n_frags * 6 + sys->n_charges * 3, sizeof(double));
 	state->ff = NULL;
 
+
 	if (cfg_get_bool(cfg, "enable_ff")) {
 		if ((state->ff = ff_create()) == NULL)
 			error("cannot create ff object");

efpmd/src/main.c

@@ -34,6 +34,7 @@ void sim_sp(struct state *state)
 
 	print_geometry(state->efp);
 	compute_energy(state, false);
+    print_pair_energy(state);
 	print_energy(state);
 
 	msg("SINGLE POINT ENERGY JOB COMPLETED SUCCESSFULLY\n");

efpmd/src/sp.c

@@ -0,0 +1,68 @@
+// !$*UTF8*$!
+{
+	archiveVersion = 1;
+	classes = {
+	};
+	objectVersion = 50;
+	objects = {
+
+/* Begin PBXGroup section */
+		D98CC31F240F6D7100E8763E = {
+			isa = PBXGroup;
+			children = (
+			);
+			sourceTree = "<group>";
+		};
+/* End PBXGroup section */
+
+/* Begin PBXProject section */
+		D98CC320240F6D7100E8763E /* Project object */ = {
+			isa = PBXProject;
+			attributes = {
+				LastUpgradeCheck = 1130;
+			};
+			buildConfigurationList = D98CC323240F6D7100E8763E /* Build configuration list for PBXProject "libefp" */;
+			compatibilityVersion = "Xcode 9.3";
+			developmentRegion = en;
+			hasScannedForEncodings = 0;
+			knownRegions = (
+				en,
+				Base,
+			);
+			mainGroup = D98CC31F240F6D7100E8763E;
+			projectDirPath = "";
+			projectRoot = "";
+			targets = (
+			);
+		};
+/* End PBXProject section */
+
+/* Begin XCBuildConfiguration section */
+		D98CC324240F6D7100E8763E /* Debug */ = {
+			isa = XCBuildConfiguration;
+			buildSettings = {
+			};
+			name = Debug;
+		};
+		D98CC325240F6D7100E8763E /* Release */ = {
+			isa = XCBuildConfiguration;
+			buildSettings = {
+			};
+			name = Release;
+		};
+/* End XCBuildConfiguration section */
+
+/* Begin XCConfigurationList section */
+		D98CC323240F6D7100E8763E /* Build configuration list for PBXProject "libefp" */ = {
+			isa = XCConfigurationList;
+			buildConfigurations = (
+				D98CC324240F6D7100E8763E /* Debug */,
+				D98CC325240F6D7100E8763E /* Release */,
+			);
+			defaultConfigurationIsVisible = 0;
+			defaultConfigurationName = Release;
+		};
+/* End XCConfigurationList section */
+	};
+	rootObject = D98CC320240F6D7100E8763E /* Project object */;
+}