Elec potential

efpmd/src/common.h

@@ -59,6 +59,7 @@ enum run_type {
 	RUN_TYPE_OPT,
 	RUN_TYPE_MD,
 	RUN_TYPE_EFIELD,
+	RUN_TYPE_ELPOT,
 	RUN_TYPE_GTEST,
 	RUN_TYPE_ETEST
 };

efpmd/src/efield.c

@@ -27,6 +27,7 @@
 #include "common.h"
 
 void sim_efield(struct state *state);
+void sim_elpot(struct state *state);
 
 static void
 print_field(size_t frag_idx, const struct efp_atom *atom, const double *field)
@@ -78,3 +79,54 @@ sim_efield(struct state *state)
 
 	msg("ELECTRIC FIELD JOB COMPLETED SUCCESSFULLY\n");
 }
+
+
+static void
+print_elpot(const struct efp_atom *atom, const double elpot)
+{
+     msg("    %s %12.8lf %12.8lf %12.8lf %12.8lf\n",
+        atom->label,
+        BOHR_RADIUS * atom->x,
+        BOHR_RADIUS * atom->y,
+        BOHR_RADIUS * atom->z,
+        elpot);
+    // msg("    ELEC POTENTIAL %12.8lf\n\n", elpot);
+}
+
+void
+sim_elpot(struct state *state)
+{
+    size_t n_frags;
+
+    msg("ELECTROSTATIC POTENTIAL JOB\n\n\n");
+
+    print_geometry(state->efp);
+    compute_energy(state, false);
+    print_energy(state);
+    check_fail(efp_get_frag_count(state->efp, &n_frags));
+
+    msg("COORDINATES IN ANGSTROMS, ELECTROSTATIC POTENTIAL IN ATOMIC UNITS\n");
+    msg("     ATOM            X            Y            Z        ELPOT \n\n");
+
+    for (size_t i = 0; i < n_frags; i++) {
+        double elpot;
+        struct efp_atom *atoms;
+        size_t n_atoms;
+
+        check_fail(efp_get_frag_atom_count(state->efp, i, &n_atoms));
+        atoms = xmalloc(n_atoms * sizeof(struct efp_atom));
+        check_fail(efp_get_frag_atoms(state->efp, i, n_atoms, atoms));
+
+        msg("ELECTROSTATIC POTENTIAL ON FRAGMENT %zu\n", i);
+
+        for (size_t j = 0; j < n_atoms; j++) {
+            check_fail(efp_get_elec_potential(state->efp, i, &atoms[j].x, &elpot));
+            print_elpot(atoms + j, elpot);
+        }
+        msg("\n");
+
+        free(atoms);
+    }
+
+    msg("ELECTROSTATIC POTENTIAL JOB COMPLETED SUCCESSFULLY\n");
+}

efpmd/src/main.c

@@ -36,6 +36,7 @@ void sim_hess(struct state *);
 void sim_opt(struct state *);
 void sim_md(struct state *);
 void sim_efield(struct state *);
+void sim_elpot(struct state *);
 void sim_gtest(struct state *);
 void sim_etest(struct state *);
 
@@ -56,6 +57,7 @@ static struct cfg *make_cfg(void)
 		"opt\n"
 		"md\n"
 		"efield\n"
+        "elpot\n"
 		"gtest\n"
         "etest\n",
 		(int []) { RUN_TYPE_SP,
@@ -64,10 +66,11 @@ static struct cfg *make_cfg(void)
 			   RUN_TYPE_OPT,
 			   RUN_TYPE_MD,
 			   RUN_TYPE_EFIELD,
+			   RUN_TYPE_ELPOT,
 			   RUN_TYPE_GTEST,
 			   RUN_TYPE_ETEST});
 
-	cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
+	cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_POINTS,
 		"xyzabc\n"
 		"points\n"
 		"rotmat\n",
@@ -160,20 +163,22 @@ static struct cfg *make_cfg(void)
 static sim_fn_t get_sim_fn(enum run_type run_type)
 {
 	switch (run_type) {
-	case RUN_TYPE_SP:
-		return sim_sp;
-	case RUN_TYPE_GRAD:
-		return sim_grad;
-	case RUN_TYPE_HESS:
-		return sim_hess;
-	case RUN_TYPE_OPT:
-		return sim_opt;
-	case RUN_TYPE_MD:
-		return sim_md;
-	case RUN_TYPE_EFIELD:
-		return sim_efield;
-	case RUN_TYPE_GTEST:
-		return sim_gtest;
+	    case RUN_TYPE_SP:
+		    return sim_sp;
+	    case RUN_TYPE_GRAD:
+		    return sim_grad;
+	    case RUN_TYPE_HESS:
+		    return sim_hess;
+	    case RUN_TYPE_OPT:
+		    return sim_opt;
+	    case RUN_TYPE_MD:
+		    return sim_md;
+	    case RUN_TYPE_EFIELD:
+		    return sim_efield;
+		case RUN_TYPE_ELPOT:
+            return sim_elpot;
+	    case RUN_TYPE_GTEST:
+		    return sim_gtest;
 		case RUN_TYPE_ETEST:
 		    return sim_etest;
 	}

src/efp.h

@@ -1062,6 +1062,20 @@ enum efp_result efp_get_frag_atoms(struct efp *efp, size_t frag_idx,
 enum efp_result efp_get_electric_field(struct efp *efp, size_t frag_idx,
     const double *xyz, double *field);
 
+/**
+ * Get electrostatic potential for a point on a fragment.
+ *
+ * @param[in] efp The efp structure.
+ * @param[in] frag_idx Index of a fragment. Must be a value between zero and
+ * the total number of fragments minus one.
+ * @param[in] xyz Coordinates of a point for which electric field should be
+ * computed.
+ * @param[out] elec_potential Electrostatic potential in atomic units.
+ * @return ::EFP_RESULT_SUCCESS on success or error code otherwise.
+ */
+enum efp_result efp_get_elec_potential(struct efp *efp, size_t frag_idx,
+        const double *xyz, double *elec_potential);
+
 /**
  * Convert rigid body torque to derivatives of energy by Euler angles.
  *

src/elec.h

@@ -93,6 +93,28 @@ quadrupole_sum(const double *quad, const vec_t *dr)
 	return sum;
 }
 
+static inline double
+octupole_sum(const double *oct, const vec_t *dr)
+{
+    /* order in which octupoles are stored */
+    enum { xxx = 0, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz };
+
+    double sum = 0.0;
+
+    sum += oct[xxx] * dr->x * dr->x * dr->x;
+    sum += oct[yyy] * dr->y * dr->y * dr->y;
+    sum += oct[zzz] * dr->z * dr->z * dr->z;
+    sum += oct[xxy] * dr->x * dr->x * dr->y * 3.0;
+    sum += oct[xxz] * dr->x * dr->x * dr->z * 3.0;
+    sum += oct[xyy] * dr->x * dr->y * dr->y * 3.0;
+    sum += oct[yyz] * dr->y * dr->y * dr->z * 3.0;
+    sum += oct[xzz] * dr->x * dr->z * dr->z * 3.0;
+    sum += oct[yzz] * dr->y * dr->z * dr->z * 3.0;
+    sum += oct[xyz] * dr->x * dr->y * dr->z * 6.0;
+
+    return sum;
+}
+
 double efp_charge_charge_energy(double, double, const vec_t *);
 double efp_charge_dipole_energy(double, const vec_t *, const vec_t *);
 double efp_charge_quadrupole_energy(double, const double *, const vec_t *);

src/electerms.c

@@ -26,28 +26,6 @@
 
 #include "elec.h"
 
-static double
-octupole_sum(const double *oct, const vec_t *dr)
-{
-	/* order in which octupoles are stored */
-	enum { xxx = 0, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz };
-
-	double sum = 0.0;
-
-	sum += oct[xxx] * dr->x * dr->x * dr->x;
-	sum += oct[yyy] * dr->y * dr->y * dr->y;
-	sum += oct[zzz] * dr->z * dr->z * dr->z;
-	sum += oct[xxy] * dr->x * dr->x * dr->y * 3.0;
-	sum += oct[xxz] * dr->x * dr->x * dr->z * 3.0;
-	sum += oct[xyy] * dr->x * dr->y * dr->y * 3.0;
-	sum += oct[yyz] * dr->y * dr->y * dr->z * 3.0;
-	sum += oct[xzz] * dr->x * dr->z * dr->z * 3.0;
-	sum += oct[yzz] * dr->y * dr->z * dr->z * 3.0;
-	sum += oct[xyz] * dr->x * dr->y * dr->z * 6.0;
-
-	return sum;
-}
-
 static double
 octupole_sum_xyz(const double *oct, const vec_t *dr, size_t axis)
 {

src/pol.c

@@ -114,6 +114,38 @@ get_multipole_field(const vec_t *xyz, const struct multipole_pt *mult_pt,
 	return field;
 }
 
+static double
+get_multipole_elec_potential(const vec_t *xyz, const struct multipole_pt *mult_pt,
+                    const struct swf *swf)
+{
+    double elpot = 0.0;
+
+    vec_t dr = {
+            xyz->x - mult_pt->x - swf->cell.x,
+            xyz->y - mult_pt->y - swf->cell.y,
+            xyz->z - mult_pt->z - swf->cell.z
+    };
+
+    double r = vec_len(&dr);
+    double r3 = r * r * r;
+    double r5 = r3 * r * r;
+    double r7 = r5 * r * r;
+
+    /* charge */
+    elpot += swf->swf * mult_pt->monopole / r;
+
+    /* dipole */
+    elpot += swf->swf * vec_dot(&mult_pt->dipole, &dr) / r3;
+
+    /* quadrupole */
+    elpot += swf->swf * quadrupole_sum(mult_pt->quadrupole, &dr) / r5;
+
+    /* octupole */
+    elpot += swf->swf * octupole_sum(mult_pt->octupole, &dr) / r7;
+
+    return elpot;
+}
+
 static vec_t
 get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx)
 {
@@ -1259,3 +1291,78 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz,
 	*((vec_t *)field) = elec_field;
 	return EFP_RESULT_SUCCESS;
 }
+
+EFP_EXPORT enum efp_result
+efp_get_elec_potential(struct efp *efp, size_t frag_idx, const double *xyz,
+                       double *elec_potential)
+{
+    assert(efp);
+    assert(frag_idx < efp->n_frag);
+    assert(xyz);
+    assert(elec_potential);
+
+    const struct frag *frag = efp->frags + frag_idx;
+    double elpot = 0.0;
+
+    for (size_t i = 0; i < efp->n_frag; i++) {
+        if (i == frag_idx || efp_skip_frag_pair(efp, i, frag_idx))
+            continue;
+
+        const struct frag *fr_i = efp->frags + i;
+        struct swf swf = efp_make_swf(efp, fr_i, frag, 0);
+        if (swf.swf == 0.0)
+            continue;
+
+        /* potential due to nuclei */
+        for (size_t j = 0; j < fr_i->n_atoms; j++) {
+            const struct efp_atom *at = fr_i->atoms + j;
+
+            vec_t dr = {
+                    xyz[0] - at->x - swf.cell.x,
+                    xyz[1] - at->y - swf.cell.y,
+                    xyz[2] - at->z - swf.cell.z
+            };
+
+            double r = vec_len(&dr);
+
+            elpot += swf.swf * at->znuc / r;
+        }
+
+        /* potential due to multipoles */
+        for (size_t j = 0; j < fr_i->n_multipole_pts; j++) {
+            const struct multipole_pt *mpt = fr_i->multipole_pts+j;
+            elpot += get_multipole_elec_potential((const vec_t *)xyz, mpt, &swf);
+         }
+
+        /* potential due to induced dipoles */
+        for (size_t j = 0; j < fr_i->n_polarizable_pts; j++) {
+            struct polarizable_pt *pt_i = fr_i->polarizable_pts + j;
+            size_t idx = fr_i->polarizable_offset + j;
+
+            vec_t dr = {
+                    xyz[0] - pt_i->x - swf.cell.x,
+                    xyz[1] - pt_i->y - swf.cell.y,
+                    xyz[2] - pt_i->z - swf.cell.z
+            };
+
+            double r = vec_len(&dr);
+            double r3 = r * r * r;
+
+            elpot += swf.swf * 0.5 * (vec_dot(&efp->indip[idx], &dr) + vec_dot(&efp->indipconj[idx], &dr)) / r3;
+        }
+    }
+
+    if (efp->opts.terms & EFP_TERM_AI_POL) {
+        /* elec potential due to nuclei from ab initio subsystem */
+        for (size_t i = 0; i < efp->n_ptc; i++) {
+            vec_t dr = vec_sub((const vec_t *)xyz, efp->ptc_xyz+i);
+
+            double r = vec_len(&dr);
+
+            elpot += efp->ptc[i] / r;
+        }
+    }
+
+    *elec_potential = elpot;
+    return EFP_RESULT_SUCCESS;
+}

tests/constraint_1.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy 0.0002280010
+coord xyzabc
 disp_damp tt
 elec_damp screen
 fraglib_path ../fraglib

tests/constraint_2.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy 0.0019051720
+coord xyzabc
 elec_damp screen
 disp_damp tt
 pol_damp tt

tests/disp_1a.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy -0.0000989033
+coord xyzabc
 terms disp
 disp_damp tt
 fraglib_path ../fraglib

tests/disp_1b.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy -0.0001007275
+coord xyzabc
 terms disp
 disp_damp overlap
 fraglib_path ../fraglib

tests/disp_1c.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy -0.0000980020
+coord xyzabc
 terms disp
 disp_damp off
 enable_pbc true

tests/disp_2a.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy -0.0014688094
+coord xyzabc
 terms disp
 disp_damp tt
 fraglib_path ../fraglib

tests/disp_2b.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy -0.0015801770
+coord xyzabc
 terms disp
 disp_damp overlap
 fraglib_path ../fraglib

tests/efield.in

@@ -1,4 +1,5 @@
 run_type efield
+coord xyzabc
 fraglib_path ../fraglib
 
 fragment h2o_l

tests/elec_1a.in

@@ -1,5 +1,6 @@
 run_type gtest
 ref_energy 0.0002900482
+coord xyzabc
 terms elec
 elec_damp screen
 fraglib_path ../fraglib