working but uncleaned revolutionized version

Author: slipchenko

efpmd/src/main.c

@@ -152,7 +152,7 @@ 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_int(cfg, "ligand", -100);
     cfg_add_bool(cfg, "enable_pairwise", false);
     cfg_add_bool(cfg, "print_pbc", false);
     cfg_add_bool(cfg, "symmetry", false);
@@ -291,7 +291,7 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 
 	if (opts.xr_cutoff == 0.0) {
 	    opts.xr_cutoff = opts.swf_cutoff;
-	    printf("xr_cutoff is set to %lf \n\n\n", opts.xr_cutoff*0.52917721092);
+	    printf("xr_cutoff is set to %lf \n\n", opts.xr_cutoff*0.52917721092);
 	}
 
 	enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord");
@@ -305,12 +305,8 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 	else
 		add_potentials(efp, cfg, sys);
 
-	for (size_t i = 0; i < sys->n_frags; i++) {
-	    if (i == opts.ligand)
-            check_fail(efp_add_ligand(efp, sys->frags[i].name));
-	    else
-            check_fail(efp_add_fragment(efp, sys->frags[i].name));
-    }
+	for (size_t i = 0; i < sys->n_frags; i++)
+	    check_fail(efp_add_fragment(efp, sys->frags[i].name));
 
 	if (sys->n_charges > 0) {
 		double q[sys->n_charges];
@@ -335,7 +331,10 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 	if (cfg_get_bool(cfg, "enable_ff"))
 		opts.terms &= ~(EFP_TERM_ELEC | EFP_TERM_POL | EFP_TERM_DISP | EFP_TERM_XR);
 
-	check_fail(efp_set_opts(efp, &opts));
+    if (opts.enable_pairwise)
+        check_fail(efp_add_ligand(efp, opts.ligand));
+
+    check_fail(efp_set_opts(efp, &opts));
 	check_fail(efp_prepare(efp));
     check_fail(efp_set_symmlist(efp));
 

efpmd/src/parse.c

@@ -253,6 +253,10 @@ struct sys *parse_input(struct cfg *cfg, const char *path)
 	if (cfg_get_bool(cfg, "enable_pbc"))
 		cfg_set_bool(cfg, "enable_cutoff", true);
 
+	// turn off pairwise calculations if ligand is not set
+	if (cfg_get_int(cfg, "ligand") == -100 && cfg_get_bool(cfg, "enable_pairwise") == true)
+	    error("Specify ligand for pairwise calculations");
+
 	check_cfg(cfg);
 	efp_stream_close(stream);
 

src/efp.c

@@ -201,11 +201,13 @@ free_ligand(struct ligand *ligand)
     if (!ligand)
         return;
 
-    free_frag(ligand->ligand_frag);
+    //free_frag(&ligand->ligand_frag);
     for (size_t i=0; i<ligand->n_ligand_pts; i++) {
-        free(ligand->ligand_pts[i].fragment_field)
+        if (ligand->ligand_pts[i].fragment_field)
+            free(ligand->ligand_pts[i].fragment_field);
     }
-    free(ligand->ligand_pts);
+    if (ligand->ligand_pts)
+        free(ligand->ligand_pts);
 }
 
 static enum efp_result
@@ -316,6 +318,7 @@ copy_frag(struct frag *dest, const struct frag *src)
 	return EFP_RESULT_SUCCESS;
 }
 
+/*
 static enum efp_result
 copy_ligand(struct ligand *dest, const struct ligand *src) {
     size_t size;
@@ -343,6 +346,7 @@ copy_ligand(struct ligand *dest, const struct ligand *src) {
         }
     }
 }
+*/
 
 // updates (shifts) parameters of fragment based on coordinates of fragment atoms
 static enum efp_result
@@ -515,11 +519,11 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to,
 
 					/* */
 					if (efp->opts.enable_pairwise) {
-                        if (i == efp->opts.ligand) {
+                        if (i == efp->ligand_index) {
                             efp->pair_energies[fr_j].exchange_repulsion = exr;
                             efp->pair_energies[fr_j].charge_penetration = ecp;
                         }
-                        if (fr_j == efp->opts.ligand) {
+                        if (fr_j == efp->ligand_index) {
                             efp->pair_energies[i].exchange_repulsion = exr;
                             efp->pair_energies[i].charge_penetration = ecp;
                         }
@@ -531,9 +535,9 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to,
 					e_elec += e_elec_tmp;
 					/* */
 					if (efp->opts.enable_pairwise) {
-                        if (i == efp->opts.ligand) 
+                        if (i == efp->ligand_index)
                             efp->pair_energies[fr_j].electrostatic = e_elec_tmp;                       
-                        if (fr_j == efp->opts.ligand) 
+                        if (fr_j == efp->ligand_index)
                             efp->pair_energies[i].electrostatic = e_elec_tmp;
                     }
 				}
@@ -1181,6 +1185,7 @@ efp_prepare(struct efp *efp)
 		efp->n_polarizable_pts += efp->frags[i].n_polarizable_pts;
 	}
 
+	/*
 	efp->n_fragment_field_pts = 0;
 	if (efp->opts.enable_pairwise && efp->opts.ligand != -1) {
 		size_t ligand_idx = efp->opts.ligand;
@@ -1190,6 +1195,7 @@ efp_prepare(struct efp *efp)
 		}
 	}
 	efp->fragment_field = (vec_t *)calloc(efp->n_fragment_field_pts, sizeof(vec_t));
+*/
 
 	efp->grad = (six_t *)calloc(efp->n_frag, sizeof(six_t));
 	efp->skiplist = (char *)calloc(efp->n_frag * efp->n_frag, 1);
@@ -1759,6 +1765,9 @@ efp_shutdown(struct efp *efp)
         free_frag(efp->lib_current[i]);
         free(efp->lib_current[i]);
     }
+
+    free_ligand(efp->ligand);
+    free(efp->ligand);
 	free(efp->frags);
 	free(efp->lib);
     free(efp->lib_current);
@@ -1769,7 +1778,7 @@ efp_shutdown(struct efp *efp)
 	free(efp->ai_orbital_energies);
 	free(efp->ai_dipole_integrals);
 	free(efp->skiplist);
-	free(efp->fragment_field);
+	// free(efp->fragment_field);
 	free(efp->pair_energies);
     free(efp->symmlist);
 	free(efp);
@@ -1872,13 +1881,45 @@ efp_add_fragment(struct efp *efp, const char *name)
 	return EFP_RESULT_SUCCESS;
 }
 
+
 EFP_EXPORT enum efp_result
-efp_add_ligand(struct efp *efp, const char *name) {
+efp_add_ligand(struct efp *efp, int ligand_index) {
+
     assert(efp);
-    assert(name);
+    assert(ligand_index >= -1);
+
+    // ligand_index=-1 means ligand is QM; -100 is default when ligand is not specified
+    // 0 and larger: ligand is a fragment with with index
+    if (ligand_index > -1) {
+        assert( (size_t)ligand_index < efp->n_frag);
+        efp->ligand_index = (size_t)ligand_index;
+
+        efp->ligand = (struct ligand *) calloc(1, sizeof(struct ligand));
+        if (efp->ligand == NULL)
+            return EFP_RESULT_NO_MEMORY;
+
+        struct ligand *lig = efp->ligand;
 
-    check_fail(efp_add_fragment(efp, name));
-    check_fail(make_ligand(efp));
+        // copy_frag(lig->ligand_frag, &efp->frags[ligand_index]);
+        lig->ligand_frag = &efp->frags[ligand_index];
+        lig->n_ligand_pts = efp->frags[ligand_index].n_polarizable_pts;
+
+        size_t size;
+        size = lig->n_ligand_pts * sizeof(struct ligand_pt);
+        lig->ligand_pts = (struct ligand_pt *) malloc(size);
+        if (lig->ligand_pts == NULL)
+            return EFP_RESULT_NO_MEMORY;
+
+        for (size_t i = 0; i < lig->n_ligand_pts; i++) {
+            struct ligand_pt *pt = lig->ligand_pts + i;
+            size = efp->n_frag * sizeof(vec_t);
+            pt->n_frag = efp->n_frag;
+            pt->fragment_field = (vec_t *) malloc(size);
+            if (!pt->fragment_field)
+                return EFP_RESULT_NO_MEMORY;
+        }
+    }
+    return EFP_RESULT_SUCCESS;
 }
 
 EFP_EXPORT enum efp_result
@@ -2303,6 +2344,17 @@ print_pol_pt(struct efp *efp, size_t frag_index, size_t pol_index) {
     printf("\n");
 }
 
+void print_ligand(struct efp *efp, size_t frag_index) {
+    if (! efp->opts.enable_pairwise)
+        return;
+    if (efp->ligand_index == frag_index) {
+        printf("Ligand index %d\n", efp->ligand_index);
+        if (efp->ligand_index > -1) {
+            printf("Number of ligand points %d", efp->ligand->n_ligand_pts);
+        }
+    }
+}
+
 void
 print_frag_info(struct efp *efp, size_t frag_index) {
     struct frag *fr = efp->frags + frag_index;
@@ -2320,5 +2372,7 @@ print_frag_info(struct efp *efp, size_t frag_index) {
         print_pol_pt(efp, frag_index, i);
     }
 
+    print_ligand(efp, frag_index);
+
     printf("\n");
 }

src/efp.h

@@ -345,6 +345,14 @@ enum efp_result efp_add_potential(struct efp *efp, const char *path);
  */
 enum efp_result efp_add_fragment(struct efp *efp, const char *name);
 
+/**
+ * Add a ligand fragment to teh system
+ * @param[in] efp
+ * @param[in] ligand_index Index of the ligand in the fragment list
+ * @return ::EFP_RESULT_SUCCESS on success or error code otherwise.
+ */
+enum efp_result efp_add_ligand(struct efp *efp, int ligand_index);
+
 /**
  * Prepare the calculation.
  *
@@ -1368,6 +1376,13 @@ void print_mult_pt(struct efp *efp, size_t frag_index, size_t pt_index);
  */
 void print_pol_pt(struct efp *efp, size_t frag_index, size_t pol_index);
 
+/**
+ * Prints information about ligand if any
+ * @param efp
+ * @param frag_index index of fragment
+ */
+void print_ligand(struct efp *efp, size_t frag_index);
+
 /**
  * Prints information on fragment
  * @param efp