Separate cutoff for exrep and turning off interactions based on swf

Author: slipchenko

efpmd/src/main.c

@@ -110,6 +110,7 @@ static struct cfg *make_cfg(void)
 	cfg_add_string(cfg, "efp_params_file", "params.efp");
 	cfg_add_bool(cfg, "enable_cutoff", false);
 	cfg_add_double(cfg, "swf_cutoff", 10.0);
+    cfg_add_double(cfg, "xr_cutoff", 0.0);
 	cfg_add_int(cfg, "max_steps", 100);
 	cfg_add_int(cfg, "multistep_steps", 1);
 	cfg_add_string(cfg, "fraglib_path", FRAGLIB_PATH);
@@ -253,11 +254,17 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 		.enable_pbc = cfg_get_bool(cfg, "enable_pbc"),
 		.enable_cutoff = cfg_get_bool(cfg, "enable_cutoff"),
 		.swf_cutoff = cfg_get_double(cfg, "swf_cutoff"),
+		.xr_cutoff = cfg_get_double(cfg, "xr_cutoff"),
         .enable_pairwise = cfg_get_bool(cfg, "enable_pairwise"), 
         .ligand = cfg_get_int(cfg, "ligand"),
         .print_pbc = cfg_get_bool(cfg, "print_pbc")
 	};
 
+	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);
+	}
+
 	enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord");
 	struct efp *efp = efp_create();
 
@@ -393,7 +400,9 @@ static void convert_units(struct cfg *cfg, struct sys *sys)
 		cfg_get_double(cfg, "pressure") * BAR_TO_AU);
 	cfg_set_double(cfg, "swf_cutoff",
 		cfg_get_double(cfg, "swf_cutoff") / BOHR_RADIUS);
-	cfg_set_double(cfg, "num_step_dist",
+    cfg_set_double(cfg, "xr_cutoff",
+            cfg_get_double(cfg, "xr_cutoff") / BOHR_RADIUS);
+    cfg_set_double(cfg, "num_step_dist",
 		cfg_get_double(cfg, "num_step_dist") / BOHR_RADIUS);
 
 	size_t n_convert = (size_t []) {

src/disp.c

@@ -280,27 +280,32 @@ efp_frag_frag_disp(struct efp *efp, size_t frag_i, size_t frag_j,
 	struct frag *fr_i = efp->frags + frag_i;
 	struct frag *fr_j = efp->frags + frag_j;
 
-	size_t n_disp_i = fr_i->n_dynamic_polarizable_pts;
-	size_t n_disp_j = fr_j->n_dynamic_polarizable_pts;
-
-	struct swf swf = efp_make_swf(efp, fr_i, fr_j);
-
-	for (size_t ii = 0, idx = 0; ii < n_disp_i; ii++)
-		for (size_t jj = 0; jj < n_disp_j; jj++, idx++)
-			energy += point_point_disp(efp, frag_i, frag_j, ii, jj,
-			    s[idx], ds[idx], &swf);
-
-	vec_t force = {
-		swf.dswf.x * energy,
-		swf.dswf.y * energy,
-		swf.dswf.z * energy
-	};
-
-	six_atomic_add_xyz(efp->grad + frag_i, &force);
-	six_atomic_sub_xyz(efp->grad + frag_j, &force);
-	efp_add_stress(&swf.dr, &force, &efp->stress);
-
-	return energy * swf.swf;
+    struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
+    // skip calculations if distance between fragments is too big...
+    if (swf.swf == 0.0) {
+        return 0.0;
+    }
+    else {
+        size_t n_disp_i = fr_i->n_dynamic_polarizable_pts;
+        size_t n_disp_j = fr_j->n_dynamic_polarizable_pts;
+
+
+        for (size_t ii = 0, idx = 0; ii < n_disp_i; ii++)
+            for (size_t jj = 0; jj < n_disp_j; jj++, idx++)
+                energy += point_point_disp(efp, frag_i, frag_j, ii, jj,
+                                           s[idx], ds[idx], &swf);
+
+        vec_t force = {
+                swf.dswf.x * energy,
+                swf.dswf.y * energy,
+                swf.dswf.z * energy
+        };
+
+        six_atomic_add_xyz(efp->grad + frag_i, &force);
+        six_atomic_sub_xyz(efp->grad + frag_j, &force);
+        efp_add_stress(&swf.dr, &force, &efp->stress);
+        return energy * swf.swf;
+    }
 }
 
 void

src/efp.c

@@ -288,6 +288,12 @@ check_opts(const struct efp_opts *opts)
 			return EFP_RESULT_FATAL;
 		}
 	}
+    if (opts->enable_cutoff) {
+        if (opts->swf_cutoff < opts->xr_cutoff) {
+            efp_log("exchange-repulsion cutoff is smaller than interaction cutoff");
+            return EFP_RESULT_FATAL;
+        }
+    }
 	return EFP_RESULT_SUCCESS;
 }
 

src/efp.h

@@ -150,7 +150,9 @@ struct efp_opts {
 	int enable_cutoff;
 	/** Cutoff distance for fragment-fragment interactions. */
 	double swf_cutoff;
-	/** Enable ligand-fragment energy decomposition from total system */ 
+    /** Cutoff distance for exchange-repulsion calculations. */
+    double xr_cutoff;
+	/** Enable ligand-fragment energy decomposition from total system */
     int enable_pairwise; 
     /** Index of ligand for enable_pairwise; default = 0; */
     int ligand;

src/elec.c

@@ -337,91 +337,97 @@ efp_frag_frag_elec(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx)
 {
 	struct frag *fr_i = efp->frags + fr_i_idx;
 	struct frag *fr_j = efp->frags + fr_j_idx;
-	struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+	struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
 	double energy = 0.0;
 
-	/* nuclei - nuclei */
-	for (size_t ii = 0; ii < fr_i->n_atoms; ii++) {
-		for (size_t jj = 0; jj < fr_j->n_atoms; jj++) {
-			struct efp_atom *at_i = fr_i->atoms + ii;
-			struct efp_atom *at_j = fr_j->atoms + jj;
-
-			vec_t dr = {
-				at_j->x - at_i->x - swf.cell.x,
-				at_j->y - at_i->y - swf.cell.y,
-				at_j->z - at_i->z - swf.cell.z
-			};
-
-			energy += efp_charge_charge_energy(at_i->znuc,
-			    at_j->znuc, &dr);
-			if (efp->do_gradient) {
-				vec_t force, add_i, add_j;
-
-				efp_charge_charge_grad(at_i->znuc, at_j->znuc,
-				    &dr, &force, &add_i, &add_j);
-				vec_scale(&force, swf.swf);
-				efp_add_force(efp->grad + fr_i_idx,
-				    CVEC(fr_i->x), CVEC(at_i->x), &force, NULL);
-				efp_sub_force(efp->grad + fr_j_idx,
-				    CVEC(fr_j->x), CVEC(at_j->x), &force, NULL);
-				efp_add_stress(&swf.dr, &force, &efp->stress);
-			}
-		}
-	}
-
-	/* nuclei - mult points */
-	for (size_t ii = 0; ii < fr_i->n_atoms; ii++) {
-		for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) {
-			energy += atom_mult_energy(efp, fr_i, fr_j,
-			    ii, jj, &swf);
-			if (efp->do_gradient) {
-				atom_mult_grad(efp, fr_i_idx, fr_j_idx,
-				    ii, jj, &swf);
-			}
-		}
-	}
-
-	/* mult points - nuclei */
-	for (size_t jj = 0; jj < fr_j->n_atoms; jj++) {
-		for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) {
-			struct swf swf2 = swf;
-
-			vec_negate(&swf2.cell);
-			vec_negate(&swf2.dr);
-			vec_negate(&swf2.dswf);
-
-			energy += atom_mult_energy(efp, fr_j, fr_i,
-			    jj, ii, &swf2);
-			if (efp->do_gradient) {
-				atom_mult_grad(efp, fr_j_idx, fr_i_idx,
-				    jj, ii, &swf2);
-			}
-		}
-	}
-
-	/* mult points - mult points */
-	for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) {
-		for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) {
-			energy += mult_mult_energy(efp, fr_i_idx, fr_j_idx,
-			    ii, jj, &swf);
-			if (efp->do_gradient) {
-				mult_mult_grad(efp, fr_i_idx, fr_j_idx,
-				    ii, jj, &swf);
-			}
-		}
-	}
-
-	vec_t force = {
-		swf.dswf.x * energy,
-		swf.dswf.y * energy,
-		swf.dswf.z * energy
-	};
-
-	six_atomic_add_xyz(efp->grad + fr_i_idx, &force);
-	six_atomic_sub_xyz(efp->grad + fr_j_idx, &force);
-	efp_add_stress(&swf.dr, &force, &efp->stress);
-
-	return energy * swf.swf;
+    // skip calculations if distance between fragments is too big...
+    if (swf.swf == 0.0) {
+        return 0.0;
+    }
+    else {
+        /* nuclei - nuclei */
+        for (size_t ii = 0; ii < fr_i->n_atoms; ii++) {
+            for (size_t jj = 0; jj < fr_j->n_atoms; jj++) {
+                struct efp_atom *at_i = fr_i->atoms + ii;
+                struct efp_atom *at_j = fr_j->atoms + jj;
+
+                vec_t dr = {
+                        at_j->x - at_i->x - swf.cell.x,
+                        at_j->y - at_i->y - swf.cell.y,
+                        at_j->z - at_i->z - swf.cell.z
+                };
+
+                energy += efp_charge_charge_energy(at_i->znuc,
+                                                   at_j->znuc, &dr);
+                if (efp->do_gradient) {
+                    vec_t force, add_i, add_j;
+
+                    efp_charge_charge_grad(at_i->znuc, at_j->znuc,
+                                           &dr, &force, &add_i, &add_j);
+                    vec_scale(&force, swf.swf);
+                    efp_add_force(efp->grad + fr_i_idx,
+                                  CVEC(fr_i->x), CVEC(at_i->x), &force, NULL);
+                    efp_sub_force(efp->grad + fr_j_idx,
+                                  CVEC(fr_j->x), CVEC(at_j->x), &force, NULL);
+                    efp_add_stress(&swf.dr, &force, &efp->stress);
+                }
+            }
+        }
+
+        /* nuclei - mult points */
+        for (size_t ii = 0; ii < fr_i->n_atoms; ii++) {
+            for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) {
+                energy += atom_mult_energy(efp, fr_i, fr_j,
+                                           ii, jj, &swf);
+                if (efp->do_gradient) {
+                    atom_mult_grad(efp, fr_i_idx, fr_j_idx,
+                                   ii, jj, &swf);
+                }
+            }
+        }
+
+        /* mult points - nuclei */
+        for (size_t jj = 0; jj < fr_j->n_atoms; jj++) {
+            for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) {
+                struct swf swf2 = swf;
+
+                vec_negate(&swf2.cell);
+                vec_negate(&swf2.dr);
+                vec_negate(&swf2.dswf);
+
+                energy += atom_mult_energy(efp, fr_j, fr_i,
+                                           jj, ii, &swf2);
+                if (efp->do_gradient) {
+                    atom_mult_grad(efp, fr_j_idx, fr_i_idx,
+                                   jj, ii, &swf2);
+                }
+            }
+        }
+
+        /* mult points - mult points */
+        for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) {
+            for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) {
+                energy += mult_mult_energy(efp, fr_i_idx, fr_j_idx,
+                                           ii, jj, &swf);
+                if (efp->do_gradient) {
+                    mult_mult_grad(efp, fr_i_idx, fr_j_idx,
+                                   ii, jj, &swf);
+                }
+            }
+        }
+
+        vec_t force = {
+                swf.dswf.x * energy,
+                swf.dswf.y * energy,
+                swf.dswf.z * energy
+        };
+
+        six_atomic_add_xyz(efp->grad + fr_i_idx, &force);
+        six_atomic_sub_xyz(efp->grad + fr_j_idx, &force);
+        efp_add_stress(&swf.dr, &force, &efp->stress);
+
+        return energy * swf.swf;
+    }
 }
 
 static void

src/pol.c

@@ -125,7 +125,9 @@ get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx)
 			continue;
 
 		const struct frag *fr_i = efp->frags + i;
-		struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+		struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
+		if (swf.swf == 0.0)
+		    continue;
 
 		/* field due to nuclei */
 		for (size_t j = 0; j < fr_i->n_atoms; j++) {
@@ -205,7 +207,9 @@ get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, size_t l
 	if (efp_skip_frag_pair(efp, ligand_idx, frag_idx))
 		return elec_field;
 
-	struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+	struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
+	if (swf.swf == 0)
+        return elec_field;
 
 	/* field due to nuclei */
 	for (size_t j = 0; j < fr_i->n_atoms; j++) {
@@ -429,7 +433,9 @@ get_induced_dipole_field(struct efp *efp, size_t frag_idx,
 			continue;
 
 		struct frag *fr_j = efp->frags + j;
-		struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+		struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
+		if (swf.swf == 0)
+		    continue;
 
 		for (size_t jj = 0; jj < fr_j->n_polarizable_pts; jj++) {
 			struct polarizable_pt *pt_j = fr_j->polarizable_pts+jj;
@@ -649,7 +655,9 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx)
 			continue;
 
 		struct frag *fr_j = efp->frags + j;
-		struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+		struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
+		if (swf.swf == 0.0)
+		    continue;
 
 		/* energy without switching applied */
 		double energy = 0.0;
@@ -917,7 +925,9 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz,
 			continue;
 
 		const struct frag *fr_i = efp->frags + i;
-		struct swf swf = efp_make_swf(efp, fr_i, frag);
+		struct swf swf = efp_make_swf(efp, fr_i, frag, 0);
+        if (swf.swf == 0.0)
+            continue;
 
 		/* field due to nuclei */
 		for (size_t j = 0; j < fr_i->n_atoms; j++) {

src/poldirect.c

@@ -65,7 +65,7 @@ get_int_mat(const struct efp *efp, size_t i, size_t j, size_t ii, size_t jj)
 	const struct frag *fr_j = efp->frags + j;
 	const struct polarizable_pt *pt_i = fr_i->polarizable_pts + ii;
 	const struct polarizable_pt *pt_j = fr_j->polarizable_pts + jj;
-	struct swf swf = efp_make_swf(efp, fr_i, fr_j);
+	struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0);
 
 	vec_t dr = {
 		pt_j->x - pt_i->x - swf.cell.x,

src/util.c

@@ -65,7 +65,7 @@ efp_skip_frag_pair(const struct efp *efp, size_t fr_i_idx, size_t fr_j_idx)
 
 struct swf
 efp_make_swf(const struct efp *efp, const struct frag *fr_i,
-    const struct frag *fr_j)
+    const struct frag *fr_j, int short_range_cutoff)
 {
 	struct swf swf;
 
@@ -97,9 +97,16 @@ efp_make_swf(const struct efp *efp, const struct frag *fr_i,
 
 	double r = vec_len(&swf.dr);
 
-	swf.swf = efp_get_swf(r, efp->opts.swf_cutoff);
-	double dswf = efp_get_dswf(r, efp->opts.swf_cutoff);
-
+	// differentiating between cutoff for exrep and other terms
+	double dswf = 0.0;
+	if (short_range_cutoff == 1) {
+        swf.swf = efp_get_swf(r, efp->opts.xr_cutoff);
+        dswf = efp_get_dswf(r, efp->opts.xr_cutoff);
+	}
+	else {
+        swf.swf = efp_get_swf(r, efp->opts.swf_cutoff);
+        dswf = efp_get_dswf(r, efp->opts.swf_cutoff);
+    }
 	swf.dswf.x = -dswf * swf.dr.x;
 	swf.dswf.y = -dswf * swf.dr.y;
 	swf.dswf.z = -dswf * swf.dr.z;

src/util.h

@@ -34,7 +34,7 @@ struct frag;
 
 int efp_skip_frag_pair(const struct efp *, size_t, size_t);
 struct swf efp_make_swf(const struct efp *, const struct frag *,
-    const struct frag *);
+    const struct frag *, int);
 int efp_check_rotation_matrix(const mat_t *);
 void efp_points_to_matrix(const double *, mat_t *);
 const struct frag *efp_find_lib(struct efp *, const char *);