Merge pull request #2 from libefp2/pbc

Pbc

Author: slipchenko

.gitignore

@@ -20,4 +20,6 @@
 /efpmd/src/efpmd
 /fraglib/params
 /bin/efpmd
+/share/libefp/fraglib/
+/include/
 tags

.idea/.gitignore

@@ -241,14 +241,25 @@ Setting `enable_pbc` to `true` also sets `enable_cutoff` to `true`.
 
 ##### Periodic Box Size
 
-`periodic_box <x> <y> <z>`
+`periodic_box <x> <y> <z> <alpha> <beta> <gamma>`
 
-Default value: `30.0 30.0 30.0`
+Default value: `30.0 30.0 30.0 90.0 90.0 90.0`
 
-Unit: Angstrom
+Unit: Angstroms, degrees
 
+If only three values are given, the angles are set to 90 degrees (orthogonal box). 
+Non-orthogonal PBC is implemented only for single-point energy calculations.
 The smallest box dimension must be greater than `2 * swf_cutoff`.
 
+##### Print PBC coordinates
+
+`print_pbc [true/false]`
+
+Default value: `false`
+
+Prints coordinates of the system contained in a single periodic cell around 
+a fragment specified by `ligand` keyword.
+
 ### Geometry optimization related parameters
 
 ##### Optimization tolerance

.idea/libefp_pbc.iml

@@ -140,6 +140,61 @@ void print_geometry(struct efp *efp)
 	msg("\n\n");
 }
 
+void print_geometry_pbc(struct efp *efp, int ligand)
+{
+    size_t n_frags;
+    check_fail(efp_get_frag_count(efp, &n_frags));
+
+    msg("    GEOMETRY IN PBC CELL (ANGSTROMS)\n\n");
+
+    double bx[6];
+
+    check_fail(efp_get_periodic_box(efp, bx));
+    six_t box = {bx[0],bx[1],bx[2],bx[3],bx[4],bx[5]};
+    //six_t box = {10.66, 12.03, 10.872, 90.0, 115.83, 90.0};
+
+    double lig_com[6];
+    check_fail(efp_get_frag_xyzabc(efp,ligand,lig_com));
+    for (size_t i = 0; i < n_frags; i++) {
+
+        vec_t cell = {0.0,0.0,0.0};
+        if (i != ligand) {
+            double frag_com[6];
+            check_fail(efp_get_frag_xyzabc(efp, i, frag_com));
+            vec_t dr = {frag_com[0] - lig_com[0], frag_com[1] - lig_com[1], frag_com[2] - lig_com[2]};
+
+            if (box.a == 90.0 && box.b == 90.0 && box.c == 90.0) {
+                cell.x = box.x * round(dr.x / box.x);
+                cell.y = box.y * round(dr.y / box.y);
+                cell.z = box.z * round(dr.z / box.z);
+                //dr = vec_sub(&dr, &cell);
+            } else {
+                cart_to_frac(box, &dr);
+                cell.x = round(dr.x);
+                cell.y = round(dr.y);
+                cell.z = round(dr.z);
+                //dr = vec_sub(&dr, &cell);
+                frac_to_cart(box, &cell);
+            }
+        }
+        size_t n_atoms;
+        check_fail(efp_get_frag_atom_count(efp, i, &n_atoms));
+
+        struct efp_atom atoms[n_atoms];
+        check_fail(efp_get_frag_atoms(efp, i, n_atoms, atoms));
+
+            for (size_t a = 0; a < n_atoms; a++) {
+            double x = (atoms[a].x - cell.x)* BOHR_RADIUS;
+            double y = (atoms[a].y - cell.y)* BOHR_RADIUS;
+            double z = (atoms[a].z - cell.z)* BOHR_RADIUS;
+
+            msg("%-16s %12.6lf %12.6lf %12.6lf\n", atoms[a].label, x, y, z);
+        }
+    }
+
+    msg("\n\n");
+}
+
 void print_energy(struct state *state)
 {
 	struct efp_energy energy;
@@ -272,104 +327,103 @@ 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;
+    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);
         }
-        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;
+        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);
         }
-        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");
-	}
+        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");
 }
 
-
-vec_t box_from_str(const char *str)
+six_t box_from_str(const char *str)
 {
-	vec_t box;
+	six_t box;
 
-	if (sscanf(str, "%lf %lf %lf", &box.x, &box.y, &box.z) < 3)
+	if (sscanf(str, "%lf %lf %lf %lf %lf %lf ", &box.x, &box.y, &box.z, &box.a, &box.b, &box.c) < 3)
 		error("incorrect box format");
 
-	vec_scale(&box, 1.0 / BOHR_RADIUS);
+	box.x *= 1.0 / BOHR_RADIUS;
+	box.y *= 1.0 / BOHR_RADIUS;
+	box.z *= 1.0 / BOHR_RADIUS;
+	// vec_scale(&box, 1.0 / BOHR_RADIUS);
 	return box;
 }
 

.idea/modules.xml

@@ -61,7 +61,6 @@ enum run_type {
 	RUN_TYPE_EFIELD,
 	RUN_TYPE_GTEST,
 };
-//RUN_TYPE_PAIRWISE
 
 enum ensemble_type {
 	ENSEMBLE_TYPE_NVE,
@@ -108,6 +107,7 @@ void *xrealloc(void *, size_t);
 
 void print_vec(const double *);
 void print_geometry(struct efp *);
+void print_geometry_pbc(struct efp *, int);
 void print_energy(struct state *);
 void print_gradient(struct state *);
 void print_fragment(const char *, const double *, const double *);
@@ -119,7 +119,7 @@ void print_pair_energy(struct state *);
 void check_fail(enum efp_result);
 void compute_energy(struct state *, bool);
 struct sys *parse_input(struct cfg *, const char *);
-vec_t box_from_str(const char *);
+six_t box_from_str(const char *);
 int efp_strcasecmp(const char *, const char *);
 int efp_strncasecmp(const char *, const char *, size_t);