nonorthogonal PBC and PBC print-out

Author: slipchenko

src/efp.c

@@ -32,6 +32,7 @@
 #include "elec.h"
 #include "private.h"
 #include "stream.h"
+#include "util.h"
 
 static void
 update_fragment(struct frag *frag)
@@ -837,35 +838,32 @@ EFP_EXPORT enum efp_result
 efp_set_periodic_box(struct efp *efp, double x, double y, double z, double alpha, double beta, double gamma)
 {
 	assert(efp);
-
-	if (x < 2.0 * efp->opts.swf_cutoff ||
-	    y < 2.0 * efp->opts.swf_cutoff ||
-	    z < 2.0 * efp->opts.swf_cutoff) {
-		efp_log("periodic box dimensions must be at least twice "
-		    "the switching function cutoff");
-		return EFP_RESULT_FATAL;
-	}
-
-	if (alpha == 0.0 && beta == 0.0 && gamma == 0.0) {
+	if (alpha < 0.01) {
         // assigning default angles of 90.0 degrees
+        //printf("\n assigning angles to 90.0 \n");
         alpha = 90.0;
         beta = 90.0;
         gamma = 90.0;
     }
 
-    if (alpha == 0.0 || beta == 0.0 || gamma == 0.0) {
-        efp_log("periodic box angle cannot be zero");
-        return EFP_RESULT_FATAL;
-	}
-
 	efp->box.x = x;
 	efp->box.y = y;
 	efp->box.z = z;
 	efp->box.a = alpha;
 	efp->box.b = beta;
 	efp->box.c = gamma;
 
-	return EFP_RESULT_SUCCESS;
+	double max_cut = max_cutoff(efp->box);
+	double cut = efp->opts.swf_cutoff;
+	if (cut > max_cut) {
+	    printf("\n Maximum allowed cutoff is %lf ", max_cut*0.52917721092);
+        printf("\n Given cutoff is %lf \n", cut*0.52917721092);
+        efp_log("periodic box dimensions must be at least twice "
+                "the switching function cutoff");
+        return EFP_RESULT_FATAL;
+	}
+
+    return EFP_RESULT_SUCCESS;
 }
 
 EFP_EXPORT enum efp_result
@@ -874,7 +872,7 @@ efp_get_periodic_box(struct efp *efp, double *xyzabc)
 	assert(efp);
 	assert(xyzabc);
 
-	xyzabc[0] = efp->box.x;
+    xyzabc[0] = efp->box.x;
 	xyzabc[1] = efp->box.y;
 	xyzabc[2] = efp->box.z;
     xyzabc[3] = efp->box.a;

src/efp.h

@@ -153,7 +153,9 @@ struct efp_opts {
 	/** Enable ligand-fragment energy decomposition from total system */ 
     int enable_pairwise; 
     /** Index of ligand for enable_pairwise; default = 0; */
-    int ligand;     
+    int ligand;
+    /** Prints fragment coordinates rearranged around ligand. Applicable for periodic simulations only. */
+    int print_pbc;
 };
 
 /** EFP energy terms. */

src/mathutil.h

@@ -382,7 +382,7 @@ matrix_to_euler(const mat_t *rotmat, double *ea, double *eb, double *ec)
 	*ec = c;
 }
 
-static inline vec_t *
+static inline void
 cart_to_frac(const six_t box, vec_t *dr)
 {
     double radian = 180.0/PI;
@@ -397,10 +397,9 @@ cart_to_frac(const six_t box, vec_t *dr)
     dr->z = (dr->z/gamma_term) / box.z;
     dr->y = ((dr->y - dr->z * box.z * beta_term)/gamma_sin) / box.y;
     dr->x = (dr->x - dr->y * box.y * gamma_cos - dr->z * box.z * beta_cos) / box.x;
-    return dr;
 }
 
-static inline vec_t *
+static inline void
 frac_to_cart(const six_t box, vec_t *dr) {
     double radian = 180.0/PI;
     double alpha_cos = cos(box.a / radian);
@@ -411,10 +410,9 @@ frac_to_cart(const six_t box, vec_t *dr) {
     double beta_term = (alpha_cos - beta_cos*gamma_cos) / gamma_sin;
     double gamma_term = sqrt(beta_sin*beta_sin - beta_term*beta_term);
 
-    dr->z = dr->z * gamma_term * box.z;
-    dr->y = dr->y * box.y * gamma_sin + dr->z * box.z * beta_term;
     dr->x = dr->x * box.x + dr->y * box.y * gamma_cos + dr->z * box.z * beta_cos;
-    return dr;
+    dr->y = dr->y * box.y * gamma_sin + dr->z * box.z * beta_term;
+    dr->z = dr->z * gamma_term * box.z;
 }
 
 #endif /* LIBEFP_MATH_UTIL_H */

src/util.c

@@ -295,3 +295,46 @@ efp_strncasecmp(const char *s1, const char *s2, size_t n)
 	}
 	return 0;
 }
+
+void
+find_plane(const vec_t pt1, const vec_t pt2, const vec_t pt3, vec_t *normal, double d) {
+    // determines plane by three points: ax + by + cz + d = 0, where a,b,c are given by a normal vector
+    vec_t vec21 = vec_sub(&pt2,&pt1);
+    vec_t vec31 = vec_sub(&pt3,&pt1);
+    *normal = vec_cross(&vec21,&vec31);
+    d = - vec_dot(normal,&pt1);
+}
+
+double
+max_cutoff(const six_t box) {
+    vec_t point000 = {0.0, 0.0, 0.0};
+    vec_t point100 = {1.0, 0.0, 0.0};
+    vec_t point010 = {0.0, 1.0, 0.0};
+    vec_t point001 = {0.0, 0.0, 1.0};
+    vec_t point = {0.5,0.5,0.5};
+    frac_to_cart(box, &point100);
+    frac_to_cart(box, &point010);
+    frac_to_cart(box, &point001);
+    frac_to_cart(box, &point);
+    // plane xy: points 000, 100, 010
+    vec_t normal;
+    double d;
+    double dist;
+    double min_dist;
+    // plane xy: points 000, 100, 010
+    find_plane(point100, point000, point010, &normal, d);
+    dist = (vec_dot(&normal, &point) + d) / vec_len(&normal);
+    dist = sqrt(dist*dist);
+    min_dist = dist;
+    // plane xz: points 000, 100, 001
+    find_plane(point100, point000, point001, &normal, d);
+    dist = (vec_dot(&normal, &point) + d) / vec_len(&normal);
+    dist = sqrt(dist*dist);
+    if (min_dist > dist) min_dist = dist;
+    // plane yz: points 000, 010, 001
+    find_plane(point010, point000, point001, &normal, d);
+    dist = (vec_dot(&normal, &point) + d) / vec_len(&normal);
+    dist = sqrt(dist*dist);
+    if (min_dist > dist) min_dist = dist;
+    return min_dist;
+}

src/util.h

@@ -48,5 +48,9 @@ void efp_rotate_t2(const mat_t *, const double *, double *);
 void efp_rotate_t3(const mat_t *, const double *, double *);
 int efp_strcasecmp(const char *, const char *);
 int efp_strncasecmp(const char *, const char *, size_t);
+// find plane by 3 points
+void find_plane(const vec_t, const vec_t, const vec_t, vec_t *, double);
+// computes maximum cut_off distance for an arbitrary periodic cell
+double max_cutoff(const six_t);
 
 #endif /* LIBEFP_UTIL_H */