match with qchem 1

Author: slipchenko

EFP2.C

@@ -1,6 +1,6 @@
 #!/bin/zsh
 
-export TORCH_SWITCH=ON
+export TORCH_SWITCH=OFF
 
 export LIBEFP_DIR="./"
 export INSTALLATION_DIR="$LIBEFP_DIR"

EFP2.h

@@ -1729,7 +1729,7 @@ efp_compute(struct efp *efp, int do_gradient)
         efp_balance_work(efp, compute_two_body_range, NULL);
 	}
 	else {  // high-symmetry crystals
-	    if ((res = compute_two_body_crystal(efp))) {
+	    if ((res = compute_two_body_crystal(efp))){
             efp_log("compute_two_body_crystal() failure");
             return res;
         }
@@ -1739,11 +1739,11 @@ efp_compute(struct efp *efp, int do_gradient)
         efp_log("efp_compute_pol() failure");
         return res;
     }
-	if ((res = efp_compute_ai_elec(efp))) {
+	if ((res = efp_compute_ai_elec(efp))){
         efp_log("efp_compute_ai_elec() failure");
         return res;
     }
-	if ((res = efp_compute_ai_disp(efp))) {
+	if ((res = efp_compute_ai_disp(efp))){
         efp_log("efp_compute_ai_disp() failure");
         return res;
     }
@@ -1928,7 +1928,8 @@ efp_get_multipole_values(struct efp *efp, double *mult)
 		for (size_t j = 0; j < frag->n_multipole_pts; j++) {
 			struct multipole_pt *pt = frag->multipole_pts + j;
 
-			*mult++ = pt->monopole;
+            // monopole and nuclear charges are added together as of 2025: LVS
+			*mult++ = pt->monopole + pt->znuc;
 
 			*mult++ = pt->dipole.x;
 			*mult++ = pt->dipole.y;
@@ -1943,6 +1944,24 @@ efp_get_multipole_values(struct efp *efp, double *mult)
 	return EFP_RESULT_SUCCESS;
 }
 
+EFP_EXPORT enum efp_result
+efp_get_mono_values(struct efp *efp, double *monopoles)
+{
+    assert(efp);
+    assert(monopoles);
+
+    for (size_t i = 0; i < efp->n_frag; i++) {
+        struct frag *frag = efp->frags + i;
+
+        for (size_t j = 0; j < frag->n_multipole_pts; j++) {
+            struct multipole_pt *pt = frag->multipole_pts + j;
+
+            *monopoles++ = pt->monopole;
+		}
+	}
+	return EFP_RESULT_SUCCESS;
+}
+
 EFP_EXPORT enum efp_result
 efp_get_dipole_values(struct efp *efp, double *dipoles)
 {
@@ -2659,6 +2678,28 @@ save_ai_field_pol_pt(struct efp *efp, struct efp_pol_pt *pol_pt, size_t frag_idx
     return EFP_RESULT_SUCCESS;
 }
 
+EFP_EXPORT enum efp_result
+save_ai_field(struct efp *efp, double *field) {
+    assert(efp);
+    assert(field);
+
+    size_t k = 0;
+    for (size_t i = 0; i < efp->n_frag; i++) {
+        struct frag *frag = efp->frags + i;
+
+        for (size_t j = 0; j < frag->n_polarizable_pts; j++) {
+            struct polarizable_pt *pt = frag->polarizable_pts + j;
+
+            pt->elec_field_wf.x = -field[k];
+            pt->elec_field_wf.y = -field[k+1];
+            pt->elec_field_wf.z = -field[k+2];
+            k = k+3;
+
+            //print_pol_pt(efp,i,j);
+        }
+    }
+    return EFP_RESULT_SUCCESS;
+}
 
 EFP_EXPORT void
 efp_torque_to_derivative(const double *euler, const double *torque,
@@ -2688,7 +2729,7 @@ efp_banner(void)
 	static const char banner[] =
 		"LIBEFP ver. " LIBEFP_VERSION_STRING "\n"
 		"Copyright (c) 2012-2017 Ilya Kaliman\n"
-        "              2018-2022 Lyudmila Slipchenko\n"
+        "              2018-2025 Lyudmila Slipchenko\n"
 		"\n"
 		"Journal References:\n"
 		"  - Kaliman and Slipchenko, JCC 2013.\n"

setup.sh

@@ -177,6 +177,8 @@ struct efp_opts {
 	enum efp_elec_damp elec_damp;
 	/** Polarization damping type (see #efp_pol_damp). */
 	enum efp_pol_damp pol_damp;
+	/** Polarization damping value (see #efp_pol_damp_tt_value). */
+	double pol_damp_tt_value;
 	/** Driver used to find polarization induced dipoles. */
 	enum efp_pol_driver pol_driver;
 	/** Enable periodic boundary conditions if nonzero. */
@@ -994,14 +996,27 @@ enum efp_result efp_get_multipole_coordinates(struct efp *efp, double *xyz);
  */
 enum efp_result efp_get_multipole_values(struct efp *efp, double *mult);
 
+/**
+ *
+ * @param efp The efp structure.
+ * @param[out] monopoles Array with all efp monopoles.
+ *
+ * The size of the \p monopoles array must be at least [\p n_mult] elements
+ * where \p n_mult is the value returned by the ::efp_get_multipole_count function.
+ * Only monopoles (no nuclear charges) are returned by this funciton.
+ *
+ * \return ::EFP_RESULT_SUCCESS on success or error code otherwise.
+ */
+enum efp_result efp_get_mono_values(struct efp *efp, double *monopoles);
+
 /**
  * Get electrostatics dipoles from EFP fragments.
  *
  * \param[in] efp The efp structure.
  *
  * \param[out] dipoles Array with all efp dipoles.
  *
- * The size of the \p mult array must be at least [3 * \p n_mult] elements
+ * The size of the \p dipoles array must be at least [3 * \p n_mult] elements
  * where \p n_mult is the value returned by the ::efp_get_multipole_count function.
  *
  * \return ::EFP_RESULT_SUCCESS on success or error code otherwise.
@@ -1402,6 +1417,14 @@ efp_get_frag_pol_pt(struct efp *efp, size_t frag_idx, size_t pt_idx,
 enum efp_result
 save_ai_field_pol_pt(struct efp *efp, struct efp_pol_pt *pol_pt, size_t frag_idx, size_t pt_idx);
 
+/**
+ * Saves ab initio field into polarizable points
+ * @param efp
+ * @param field pointer to the ab initio field
+ * @return
+ */
+enum efp_result save_ai_field(struct efp *efp, double *field);
+
 /**
  * Get electric field for a point on a fragment.
  *

src/efp.c

@@ -1565,7 +1565,8 @@ parse_file(struct efp *efp, struct stream *stream)
 		efp->lib[efp->n_lib - 1] = frag;
 
 		/* default value */
-		frag->pol_damp = 0.6;
+		//frag->pol_damp = 0.6; //commented out for ticket 3416
+		frag->pol_damp = efp->opts.pol_damp_tt_value;
 
 		efp_stream_next_line(stream);
 		efp_stream_next_line(stream);

src/efp.h

@@ -959,6 +959,7 @@ efp_compute_pol_energy(struct efp *efp, double *energy)
 
 	assert(energy);
 
+    printf("\nin efp_compute_pol_energy\n");
 	// counter to know when to zero out induced dipoles and static field
 	// need to be explored further
 	static int counter = 0;
@@ -1290,6 +1291,8 @@ efp_compute_pol(struct efp *efp)
 	    !(efp->opts.terms & EFP_TERM_AI_POL))
 		return EFP_RESULT_SUCCESS;
 
+    printf("\nin efp_compute_pol\n");
+
 	// this is standard non-symmetric case
 	if (! efp->opts.symmetry) {
         if ((res = efp_compute_pol_energy(efp, &efp->energy.polarization)))

src/parse.c

@@ -0,0 +1,12 @@
+run_type grad
+coord xyzabc
+terms elec pol
+elec_damp screen
+fraglib_path ../fraglib
+print 2
+
+fragment h2o_l
+   0.0   0.0   0.0   1.0   2.0   3.0
+
+fragment nh3_l
+   5.0   0.0   0.0   5.0   2.0   8.0