diff --git a/.DS_Store b/.DS_Store index 6d4ee30..6549802 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/Doxyfile b/Doxyfile index daaa01f..bf369ed 100644 --- a/Doxyfile +++ b/Doxyfile @@ -58,7 +58,7 @@ PROJECT_LOGO = # entered, it will be relative to the location where doxygen was started. If # left blank the current directory will be used. -OUTPUT_DIRECTORY = +OUTPUT_DIRECTORY = ../WEB # If the CREATE_SUBDIRS tag is set to YES then doxygen will create 4096 sub- # directories (in 2 levels) under the output directory of each output format and @@ -790,7 +790,7 @@ WARN_LOGFILE = # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING # Note: If this tag is empty the current directory is searched. -INPUT = src/efp.h +INPUT = src efpmd/src # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses diff --git a/src/efp.c b/src/efp.c index 37aa2a5..f2f56b9 100644 --- a/src/efp.c +++ b/src/efp.c @@ -502,6 +502,9 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, (void)data; + // ligand is a fragment (not QM) + bool if_pairwise = efp->opts.enable_pairwise && efp->opts.ligand > -1; + #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) reduction(+:e_elec,e_disp,e_xr,e_cp) #endif @@ -533,7 +536,7 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, e_cp += ecp; /* */ - if (efp->opts.enable_pairwise) { + if (if_pairwise) { if (i == efp->ligand_index) { efp->pair_energies[fr_j].exchange_repulsion = exr; efp->pair_energies[fr_j].charge_penetration = ecp; @@ -547,11 +550,10 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, } if (do_elec(&efp->opts) && efp->frags[i].n_multipole_pts > 0 && efp->frags[fr_j].n_multipole_pts > 0) { - e_elec_tmp = efp_frag_frag_elec(efp, - i, fr_j); + e_elec_tmp = efp_frag_frag_elec(efp, i, fr_j); e_elec += e_elec_tmp; /* */ - if (efp->opts.enable_pairwise) { + if (if_pairwise) { if (i == efp->ligand_index) efp->pair_energies[fr_j].electrostatic = e_elec_tmp; if (fr_j == efp->ligand_index) @@ -564,7 +566,7 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, i, fr_j, s, ds); e_disp += e_disp_tmp; /* */ - if (efp->opts.enable_pairwise) { + if (if_pairwise) { if (i == efp->opts.ligand) efp->pair_energies[fr_j].dispersion = e_disp_tmp; if (fr_j == efp->opts.ligand) @@ -582,6 +584,12 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, efp->energy.dispersion += e_disp; efp->energy.exchange_repulsion += e_xr; efp->energy.charge_penetration += e_cp; + + if (efp->opts.print > 0) { + printf(" In compute_two_body_range() \n"); + print_ene(&efp->energy); + print_energies(efp); + } } EFP_EXPORT enum efp_result @@ -1208,18 +1216,6 @@ efp_prepare(struct efp *efp) efp->n_polarizable_pts += efp->frags[i].n_polarizable_pts; } - /* - efp->n_fragment_field_pts = 0; - if (efp->opts.enable_pairwise && efp->opts.ligand != -1) { - size_t ligand_idx = efp->opts.ligand; - for (size_t i = 0; i < efp->n_frag; i++) { - efp->frags[i].fragment_field_offset = efp->n_fragment_field_pts; - efp->n_fragment_field_pts += efp->frags[ligand_idx].n_polarizable_pts; - } - } - efp->fragment_field = (vec_t *)calloc(efp->n_fragment_field_pts, sizeof(vec_t)); -*/ - efp->grad = (six_t *)calloc(efp->n_frag, sizeof(six_t)); efp->skiplist = (char *)calloc(efp->n_frag * efp->n_frag, 1); efp->pair_energies = (struct efp_energy *)calloc(efp->n_frag, sizeof(struct efp_energy)); @@ -1421,38 +1417,28 @@ efp_get_frag_multipole_count(struct efp *efp, size_t frag_idx, size_t *n_mult) /* EFP_EXPORT enum efp_result -efp_get_frag_mult_rank(struct efp *efp, size_t frag_idx, size_t mult_idx, size_t *rank) +efp_get_frag_rank(struct efp *efp, size_t frag_idx, size_t *rank) { assert(efp); assert(rank); assert(frag_idx < efp->n_frag); - assert(mult_idx < efp->frags[frag_idx].n_multipole_pts); struct frag *frag = efp->frags + frag_idx; - struct multipole_pt *pt = frag->multipole_pts + mult_idx; - - *rank = -1; - - for (size_t t = 0; t < 10; t++) - if (pt->octupole[t] != 0) { - *rank = 3; - return EFP_RESULT_SUCCESS; - } - for (size_t t = 0; t < 6; t++) - if (pt->quadrupole[t] != 0) { - *rank = 2; - return EFP_RESULT_SUCCESS; - } - - if (pt->dipole.x != 0 || pt->dipole.y != 0 || pt->dipole.z != 0) { - *rank = 1; - return EFP_RESULT_SUCCESS; - } - - if (pt->monopole != 0) { - *rank = 0; - return EFP_RESULT_SUCCESS; + *rank = 0; + for (size_t i=0; in_multipole_pts; i++) { + struct multipole_pt *pt = frag->multipole_pts + i; + size_t rank_tmp = 0; + if (pt->if_dip) + rank_tmp = 1; + if (pt->if_quad) + rank_tmp = 2; + if (pt->if_oct) + rank_tmp = 3; + if (rank_tmp > *rank) + *rank = rank_tmp; + if (*rank == 3) + break; } return EFP_RESULT_SUCCESS; @@ -1678,7 +1664,6 @@ efp_get_induced_dipole_conj_values(struct efp *efp, double *dip) 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; @@ -1836,7 +1821,6 @@ efp_shutdown(struct efp *efp) free(efp->ai_orbital_energies); free(efp->ai_dipole_integrals); free(efp->skiplist); - // free(efp->fragment_field); free(efp->pair_energies); free(efp->symmlist); free(efp); @@ -2007,27 +1991,6 @@ efp_create(void) return efp; } -EFP_EXPORT enum efp_result -efp_set_electron_density_field_fn(struct efp *efp, - efp_electron_density_field_fn fn) -{ - assert(efp); - - efp->get_electron_density_field = fn; - - return EFP_RESULT_SUCCESS; -} - -EFP_EXPORT enum efp_result -efp_set_electron_density_field_user_data(struct efp *efp, void *user_data) -{ - assert(efp); - - efp->get_electron_density_field_user_data = user_data; - - return EFP_RESULT_SUCCESS; -} - EFP_EXPORT enum efp_result efp_get_frag_count(struct efp *efp, size_t *n_frag) { @@ -2182,6 +2145,65 @@ efp_get_frag_mult_pt(struct efp *efp, size_t frag_idx, size_t pt_idx, return EFP_RESULT_SUCCESS; } +EFP_EXPORT enum efp_result +efp_get_frag_pol_pt(struct efp *efp, size_t frag_idx, size_t pt_idx, + struct efp_pol_pt *pol_pt) +{ + assert(efp); + assert(pol_pt); + assert(frag_idx < efp->n_frag); + assert(pt_idx < efp->frags[frag_idx].n_polarizable_pts); + + struct frag *frag; + frag = efp->frags + frag_idx; + struct polarizable_pt *pt; + pt = frag->polarizable_pts + pt_idx; + + pol_pt->x = pt->x; + pol_pt->y = pt->y; + pol_pt->z = pt->z; + pol_pt->indip[0] = pt->indip.x; + pol_pt->indip[1] = pt->indip.y; + pol_pt->indip[2] = pt->indip.z; + pol_pt->indipconj[0] = pt->indipconj.x; + pol_pt->indipconj[1] = pt->indipconj.y; + pol_pt->indipconj[2] = pt->indipconj.z; + pol_pt->indip_gs[0] = pt->indip_gs.x; + pol_pt->indip_gs[1] = pt->indip_gs.y; + pol_pt->indip_gs[2] = pt->indip_gs.z; + pol_pt->indipconj_gs[0] = pt->indipconj_gs.x; + pol_pt->indipconj_gs[1] = pt->indipconj_gs.y; + pol_pt->indipconj_gs[2] = pt->indipconj_gs.z; + pol_pt->ai_field[0] = pt->elec_field_wf.x; + pol_pt->ai_field[1] = pt->elec_field_wf.y; + pol_pt->ai_field[2] = pt->elec_field_wf.z; + + return EFP_RESULT_SUCCESS; +} + +EFP_EXPORT 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) { + assert(efp); + assert(pol_pt); + assert(frag_idx < efp->n_frag); + assert(pt_idx < efp->frags[frag_idx].n_polarizable_pts); + + struct frag *frag; + frag = efp->frags + frag_idx; + struct polarizable_pt *pt; + pt = frag->polarizable_pts + pt_idx; + + pt->elec_field_wf.x = pol_pt->ai_field[0]; + pt->elec_field_wf.y = pol_pt->ai_field[1]; + pt->elec_field_wf.z = pol_pt->ai_field[2]; + + if (efp->opts.print > 1) + print_pol_pt(efp,frag_idx,pt_idx); + + return EFP_RESULT_SUCCESS; +} + + EFP_EXPORT void efp_torque_to_derivative(const double *euler, const double *torque, double *deriv) @@ -2437,6 +2459,8 @@ print_pol_pt(struct efp *efp, size_t frag_index, size_t pol_index) { printf(" conjugated induced dipole %lf %lf %lf\n", pt->indipconj.x, pt->indipconj.y, pt->indipconj.z); printf(" old induced dipole %lf %lf %lf\n", pt->indip_old.x, pt->indip_old.y, pt->indip_old.z); printf(" old conjugated induced dipole %lf %lf %lf\n", pt->indipconj_old.x, pt->indipconj_old.y, pt->indipconj_old.z); + printf(" gr state induced dipole %lf %lf %lf\n", pt->indip_gs.x, pt->indip_gs.y, pt->indip_gs.z); + printf(" gr state conjug induced dipole %lf %lf %lf\n", pt->indipconj_gs.x, pt->indipconj_gs.y, pt->indipconj_gs.z); printf("\n"); } @@ -2470,17 +2494,6 @@ print_frag_info(struct efp *efp, size_t frag_index) { } print_ligand(efp, frag_index); - - for (int i=0; i < fr->n_multipole_pts; i++) { - struct efp_mult_pt *pt; - pt = (struct efp_mult_pt *)malloc(sizeof(struct efp_mult_pt)); - - efp_get_frag_mult_pt(efp, frag_index, i, pt); - print_efp_mult_pt(pt); - - free(pt); - } - printf("\n"); } @@ -2498,3 +2511,46 @@ print_efp_mult_pt(struct efp_mult_pt *pt) { printf(" screen0 = %lf, if_screen0 = %s\n", pt->screen0, pt->if_screen ? "true" : "false"); printf("\n"); } + +void +print_efp_pol_pt(struct efp_pol_pt *pt) { + printf("\n Polarizability point \n"); + printf(" Coordinates %lf %lf %lf\n", pt->x, pt->y, pt->z); + printf(" induced dipole %lf %lf %lf\n", pt->indip[0], pt->indip[1], pt->indip[2]); + printf(" conjug induced dipole %lf %lf %lf\n", pt->indipconj[0], pt->indipconj[1], pt->indipconj[2]); + printf(" gr state induced dipole %lf %lf %lf\n", pt->indip_gs[0], pt->indip_gs[1], pt->indip_gs[2]); + printf(" gr state conjug induced dipole %lf %lf %lf\n", pt->indipconj_gs[0], pt->indipconj_gs[1], pt->indipconj_gs[2]); + printf(" AI field %lf %lf %lf\n", pt->ai_field[0], pt->ai_field[1], pt->ai_field[2]); + printf("\n"); +} + +void print_ene(struct efp_energy *energy) { + + printf(" EFP ENERGY COMPONENTS \n"); + + printf(" ELECTROSTATIC ENERGY %lf \n", energy->electrostatic); + printf(" AI ELECTROSTATIC ENERGY %lf \n", energy->ai_electrostatic); + printf(" CHARGE PENETRATION ENERGY %lf \n", energy->charge_penetration); + printf(" POINT CHARGES ENERGY %lf \n", energy->electrostatic_point_charges); + printf(" POLARIZATION ENERGY %lf \n", energy->polarization); + printf(" EXC STATE POLARIZATION ENERGY %lf \n", energy->exs_polarization); + printf(" AI POLARIZATION ENERGY %lf \n", energy->ai_polarization); + printf(" DISPERSION ENERGY %lf \n", energy->dispersion); + printf(" AI DISPERSION ENERGY %lf \n", energy->ai_dispersion); + printf(" EXCHANGE-REPULSION ENERGY %lf \n", energy->exchange_repulsion); + double ene_sum = energy->electrostatic + energy->charge_penetration + + energy->electrostatic_point_charges + energy->polarization + + energy->dispersion + energy->exchange_repulsion; + printf(" SUM ENERGY %lf \n", ene_sum); + printf(" TOTAL ENERGY %lf \n", energy->total); + printf("\n"); +} + +void print_energies(struct efp *efp) { + printf(" --- PAIRWISE ENERGIES --- \n"); + for (size_t i=0; in_frag; i++) { + printf(" PAIR ENERGY on FRAGMENT %d %s \n", i, efp->frags[i].name); + print_ene(&efp->pair_energies[i]); + } + printf("\n"); +} diff --git a/src/efp.h b/src/efp.h index 0171c8a..e041490 100644 --- a/src/efp.h +++ b/src/efp.h @@ -41,7 +41,7 @@ extern "C" { #endif /** Version string. */ -#define LIBEFP_VERSION_STRING "1.7.2" +#define LIBEFP_VERSION_STRING "1.8.0" /** Result of an operation. */ enum efp_result { @@ -250,6 +250,18 @@ struct efp_mult_pt { bool if_screen; /**< If screen0 parameter exists and meaningful */ }; +/** Polarizability point for working with external programs */ +struct efp_pol_pt { + double x; /**< X coordinate */ + double y; /**< Y coordinate */ + double z; /**< Z coordinate */ + double indip[3]; /**< induced dipole */ + double indipconj[3]; /**< conjugated induced dipole */ + double indip_gs[3]; /**< induced dipole of the ground electronic state*/ + double indipconj_gs[3]; /**< conjugated induced dipole of the ground electronic state */ + double ai_field[3]; /** < field due to QM wavefunction */ +}; + /** * Callback function which is called by libefp to obtain electric field in the * specified points. @@ -272,8 +284,8 @@ struct efp_mult_pt { * \return The implemented function should return ::EFP_RESULT_FATAL on error * and ::EFP_RESULT_SUCCESS if the calculation has succeeded. */ -typedef enum efp_result (*efp_electron_density_field_fn)(size_t n_pt, - const double *xyz, double *field, void *user_data); +//typedef enum efp_result (*efp_electron_density_field_fn)(size_t n_pt, +// const double *xyz, double *field, void *user_data); /** * Get a human readable banner string with information about the library. @@ -408,8 +420,8 @@ enum efp_result efp_skip_fragments(struct efp *efp, size_t i, size_t j, * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -enum efp_result efp_set_electron_density_field_fn(struct efp *efp, - efp_electron_density_field_fn fn); +// enum efp_result efp_set_electron_density_field_fn(struct efp *efp, +// efp_electron_density_field_fn fn); /** * Set user data to be passed to ::efp_electron_density_field_fn. @@ -421,8 +433,8 @@ enum efp_result efp_set_electron_density_field_fn(struct efp *efp, * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -enum efp_result efp_set_electron_density_field_user_data(struct efp *efp, - void *user_data); +// enum efp_result efp_set_electron_density_field_user_data(struct efp *efp, +// void *user_data); /** * Setup arbitrary point charges interacting with EFP subsystem. @@ -829,15 +841,12 @@ enum efp_result efp_get_frag_multipole_coord(struct efp *efp, size_t frag_idx, size_t *n_mult); /** - * - * @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] mult_idx Index of a multipole point. Must be a value between 0 and - * the total number of multipoles in thios fragment minus one. - * @param[out] rank Highest rank of multipole (0 - charge, 1 - dipole, 2 - quad, 3 - oct). - * If all values are zero, rank = -1. - * @return + * Computes multipole rank of a fragment + * @param efp + * @param[in] frag_idx fragment index + * @param[out] rank Highest rank of multipoles in the fragment + * (0 - charge, 1 - dipole, 2 - quad, 3 - oct) + * @return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ // enum efp_result efp_get_frag_mult_rank(struct efp *efp, size_t frag_idx, size_t mult_idx, size_t *rank); @@ -1251,6 +1260,29 @@ enum efp_result efp_get_frag_mult_pt(struct efp *efp, size_t frag_idx, size_t pt_idx, struct efp_mult_pt *mult_pt); +/** + * + * @param efp + * @param[in] frag_idx Fragment index + * @param[in] pt_idx Polarizability point index + * @param[out] pol_pt efp_pol_pt to be returned + * @return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result +efp_get_frag_pol_pt(struct efp *efp, size_t frag_idx, size_t pt_idx, + struct efp_pol_pt *pol_pt); + +/** + * Saves ab initio field info from an efp_pol_pt structure into polarizable_pt on a fragment frag_idx + * @param efp + * @param pol_pt Polarization point + * @param frag_idx + * @param pt_idx + * @return + */ +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); + /** * Get electric field for a point on a fragment. * @@ -1441,6 +1473,24 @@ void print_frag_info(struct efp *efp, size_t frag_index); */ void print_efp_mult_pt(struct efp_mult_pt *pt); +/** + * Prints information about efp_pol_pt object + * @param pt + */ +void print_efp_pol_pt(struct efp_pol_pt *pt); + +/** + * Prints all information of efp_energy structure + * @param energy + */ +void print_ene(struct efp_energy *energy); + +/** + * Prints all information of efp->pair_energies + * @param efp + */ +void print_energies(struct efp *efp); + #ifdef __cplusplus } /* extern "C" */ #endif diff --git a/src/parse.c b/src/parse.c index 3055e01..3040cb3 100644 --- a/src/parse.c +++ b/src/parse.c @@ -46,6 +46,10 @@ static void init_multipole_pt(struct multipole_pt *pt) { pt->if_scr0 = false; } +static void init_pol_pt(struct polarizable_pt *pt) { + memset(pt, 0, sizeof(*pt)); +} + static int tok(struct stream *stream, const char *id) { @@ -464,6 +468,8 @@ parse_polarizable_pts(struct frag *frag, struct stream *stream) struct polarizable_pt *pt = frag->polarizable_pts + frag->n_polarizable_pts - 1; + // zero out all entries + init_pol_pt(pt); if (!efp_stream_advance(stream, 4)){ efp_log("parse_polarizable_pts() failure"); diff --git a/src/pol.c b/src/pol.c index af6ad6e..98a56cc 100644 --- a/src/pol.c +++ b/src/pol.c @@ -38,12 +38,6 @@ double efp_get_pol_damp_tt(double, double, double); enum efp_result efp_compute_id_direct(struct efp *); -struct id_work_data { - double conv; - vec_t *id_new; - vec_t *id_conj_new; -}; - double efp_get_pol_damp_tt(double r, double pa, double pb) { @@ -279,49 +273,6 @@ get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, int liga return elec_field; } -static enum efp_result -add_electron_density_field(struct efp *efp) -{ - enum efp_result res; - vec_t *xyz, *field; - - if (efp->get_electron_density_field == NULL) - return EFP_RESULT_SUCCESS; - - xyz = (vec_t *)malloc(efp->n_polarizable_pts * sizeof(vec_t)); - field = (vec_t *)malloc(efp->n_polarizable_pts * sizeof(vec_t)); - - for (size_t i = 0, idx = 0; i < efp->n_frag; i++) { - struct frag *frag = efp->frags + i; - - for (size_t j = 0; j < frag->n_polarizable_pts; j++, idx++) { - struct polarizable_pt *pt = frag->polarizable_pts + j; - - xyz[idx].x = pt->x; - xyz[idx].y = pt->y; - xyz[idx].z = pt->z; - } - } - - if ((res = efp->get_electron_density_field(efp->n_polarizable_pts, - (const double *)xyz, (double *)field, - efp->get_electron_density_field_user_data))) - goto error; - - for (size_t i = 0, idx = 0; i < efp->n_frag; i++) { - struct frag *frag = efp->frags + i; - - for (size_t j = 0; j < frag->n_polarizable_pts; j++, idx++) { - struct polarizable_pt *pt = frag->polarizable_pts + j; - pt->elec_field_wf = field[idx]; - } - } -error: - free(xyz); - free(field); - return res; -} - static void compute_elec_field_range(struct efp *efp, size_t from, size_t to, void *data) { @@ -381,39 +332,23 @@ compute_fragment_field_range(struct efp *efp, size_t from, size_t to, void *data } static enum efp_result -compute_elec_field(struct efp *efp) -{ - int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand != -1) ? 1 : 0; +compute_elec_field(struct efp *efp) { + enum efp_result res; efp_balance_work(efp, compute_elec_field_range, NULL); - // WARNING!!! Should it be if(efp->opts.enable_pairwise) instead ??? - if (do_pairwise) { + if (efp->opts.enable_pairwise) { // this is field due to ligand on a fragment point + // for QM ligand this is a field due to QM nuclei efp_balance_work(efp, compute_ligand_field_range, NULL); + // no contribution if ligand is QM if (efp->opts.ligand != -1) { // this is field due to fragment(s) on ligand points efp_balance_work(efp, compute_fragment_field_range, NULL); } } - // IS THIS OPENMP NEEDED HERE??? -#ifdef _OPENMP -#pragma omp parallel for schedule(dynamic) -#endif - 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++) { - frag->polarizable_pts[j].elec_field_wf = vec_zero; - } - } - - if (efp->opts.terms & EFP_TERM_AI_POL) { - if ((res = add_electron_density_field(efp))) - return res; - } - return EFP_RESULT_SUCCESS; } @@ -553,6 +488,25 @@ zero_ind_dipoles(struct efp *efp, size_t from, size_t to, void *data) } } +static void +copy_indip_gs(struct efp *efp, size_t from, size_t to, void *data) +{ + +#ifdef _OPENMP +#pragma omp parallel for schedule(dynamic) +#endif + for (size_t i = from; i < to; 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->indip_gs = pt->indip; + pt->indipconj_gs = pt->indipconj; + } + } +} + static void zero_static_field(struct efp *efp, size_t from, size_t to, void *data) { @@ -725,12 +679,12 @@ static void compute_energy_range(struct efp *efp, size_t from, size_t to, void *data) { double energy = 0.0; - // double energy1 = 0.0; - // double energy2 = 0.0; - int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand > -1) ? 1 : 0; + + // if_pairwise tells whether pairwise calculations with EFP (not QM) ligand + bool if_pairwise = efp->opts.enable_pairwise && efp->opts.ligand > -1; struct ligand *ligand; - if (do_pairwise) + if (if_pairwise) ligand = efp->ligand; #ifdef _OPENMP @@ -738,47 +692,50 @@ compute_energy_range(struct efp *efp, size_t from, size_t to, void *data) #endif for (size_t i = from; i < to; i++) { struct frag *frag = efp->frags + i; - double ene1 = 0.0, ene2 = 0.0; - // zeroing out polarization pair energies to avoid possible double counting for excited states + + // zeroing out polarization pair energies is a must efp->pair_energies[i].polarization = 0.0; + for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; energy += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf) - 0.5 * vec_dot(&pt->indip, &pt->elec_field); - // energy1 += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf); - // energy2 += -0.5 * vec_dot(&pt->indip, &pt->elec_field); - - // !!!! WARNING !!! DOES IT WORK FOR QM LIGAND ???? - if (efp->opts.enable_pairwise && i != efp->opts.ligand) { - efp->pair_energies[i].polarization += - - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); - ene2 += - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); - } + // This part is for non-QM ligand only + // Interaction of indip on this fragment with ligand's field on this fragment + if (if_pairwise) + if (i != efp->ligand_index) + efp->pair_energies[i].polarization += - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); // ligand is a QM region + // all fragments get contributions due to interaction with QM wavefunction and QM nuclei + // (ligand field contains field due to QM nuclei in the case of QM ligand) if (efp->opts.enable_pairwise && efp->opts.ligand == -1) { - efp->pair_energies[i].polarization += - 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf); - // ene1 += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf); + efp->pair_energies[i].polarization += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf); + efp->pair_energies[i].polarization += -0.5 * vec_dot(&pt->indip, &pt->ligand_field); } } - // !!!! WARNING !!! DOES IT WORK FOR QM LIGAND ???? - if (do_pairwise && i != efp->opts.ligand) { - for (size_t lp = 0; lp < ligand->n_ligand_pts; lp++) { - struct ligand_pt *lpt = ligand->ligand_pts + lp; - struct polarizable_pt *pt = ligand->ligand_frag->polarizable_pts + lp; + // this part is for non-QM ligand only + // interaction of ligand indip with the field due to this fragment + if (if_pairwise) + if ( i != efp->ligand_index ) + for (size_t lp = 0; lp < ligand->n_ligand_pts; lp++) { + struct ligand_pt *lpt = ligand->ligand_pts + lp; + struct polarizable_pt *pt = ligand->ligand_frag->polarizable_pts + lp; - efp->pair_energies[i].polarization += - - 0.5 * vec_dot(&pt->indip, &lpt->fragment_field[i]); - } - } - //printf("\n frag %d, pol pair energy1 = %lf", i, efp->pair_energies[i].polarization); + efp->pair_energies[i].polarization += + - 0.5 * vec_dot(&pt->indip, &lpt->fragment_field[i]); + } + } + + if (efp->opts.print > 1 && efp->opts.enable_pairwise) { + printf("\n Pairwise analysis from compute_energy_range() follows \n"); + print_energies(efp); } - *(double *)data += energy; + *(double *)data += energy; } static void @@ -786,13 +743,15 @@ compute_energy_correction_range(struct efp *efp, size_t from, size_t to, void *d { double energy = 0.0; + // tells whether pairwise analysis with non-QM ligand + bool if_pairwise = efp->opts.enable_pairwise && efp->opts.ligand > -1; + #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) reduction(+:energy) #endif - for (size_t i = from; i < to; i++) { struct frag *frag = efp->frags + i; - // zeroing out polarization pair energies to avoid possible double counting for excited states + // zeroing out polarization pair energies to avoid double-counting for excited states efp->pair_energies[i].exs_polarization = 0.0; for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; @@ -800,27 +759,31 @@ compute_energy_correction_range(struct efp *efp, size_t from, size_t to, void *d energy += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf) - 0.5 * vec_dot(&pt->indip, &pt->elec_field); - vec_t vec1 = vec_sub(&pt->indipconj,&pt->indipconj_old); - vec_t vec2 = vec_sub(&pt->indip, &pt->indip_old); + vec_t vec1 = vec_sub(&pt->indipconj,&pt->indipconj_gs); + vec_t vec2 = vec_sub(&pt->indip, &pt->indip_gs); vec_t ddip = vec_add(&vec1, &vec2); energy -= 0.5 * vec_dot( &ddip, &pt->elec_field_wf); - if (efp->opts.enable_pairwise && i != efp->opts.ligand) { - efp->pair_energies[i].exs_polarization += - - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); - } + // non-QM ligand. Is this part ever used? + if (if_pairwise) + if (i != efp->ligand_index) + efp->pair_energies[i].exs_polarization += + - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); // ligand is a QM region + // contributions due to excited state wavefunction and QM nuclei if (efp->opts.enable_pairwise && efp->opts.ligand == -1) { efp->pair_energies[i].exs_polarization += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf); efp->pair_energies[i].exs_polarization -= 0.5 * vec_dot(&ddip, &pt->elec_field_wf); + efp->pair_energies[i].exs_polarization += + - 0.5 * vec_dot(&pt->indip, &pt->ligand_field); } } + // subtract ground state polarization contributions efp->pair_energies[i].exs_polarization -= efp->pair_energies[i].polarization; - //printf("\n frag %d, pol pair = %lf, exc pair = %lf", i, efp->pair_energies[i].polarization, efp->pair_energies[i].exs_polarization); } *(double *)data += energy; @@ -907,22 +870,11 @@ efp_compute_pol_energy(struct efp *efp, double *energy) assert(energy); // counter to know when to zero out induced dipoles and static field + // need to be explored further static int counter = 0; - // return true (to zero out induced dipoles) only when efp_compute_pol_energy is called - // the very first time - // do not zero induced dipoles in md, optimization, qm scf iterations etc... - if ( if_clean_indip(efp, counter) ) { - efp_balance_work(efp, zero_ind_dipoles, NULL); - } - - // return true (to zero out static fields - // the very first time // think how to skip recomputing static field in qm scf iterations // check on efp->do_gradient breaks gtests... - if ( if_clean_field(efp, counter) ) - efp_balance_work(efp, zero_static_field, NULL); - if ((res = compute_elec_field(efp))) return res; @@ -951,6 +903,8 @@ efp_compute_pol_energy(struct efp *efp, double *energy) efp_balance_work(efp, compute_energy_range, energy); efp_allreduce(energy, 1); + efp_balance_work(efp, copy_indip_gs, NULL); + counter++; return EFP_RESULT_SUCCESS; @@ -963,12 +917,9 @@ efp_compute_pol_correction(struct efp *efp, double *energy) assert(energy); - //efp->energy.gs_polarization = efp->energy.polarization; - //memcpy(efp->indip_old, efp->indip, efp->n_polarizable_pts * sizeof(vec_t)); - //memcpy(efp->indipconj_old, efp->indipconj, efp->n_polarizable_pts * sizeof(vec_t)); - - if ((res = compute_elec_field(efp))) - return res; + // do not need to recompute static field for excited state correction + //if ((res = compute_elec_field(efp))) + // return res; switch (efp->opts.pol_driver) { case EFP_POL_DRIVER_ITERATIVE: @@ -1001,20 +952,8 @@ efp_compute_pol_energy_crystal(struct efp *efp, double *energy) // counter to know when to zero out induced dipoles and static field static int counter_crystal = 0; - // return true (to zero out induced dipoles) only when efp_compute_pol_energy is called - // the very first time - // do not zero induced dipoles in md, optimization, qm scf iterations etc... - if ( if_clean_indip(efp, counter_crystal) ) { - efp_balance_work(efp, zero_ind_dipoles, NULL); - } - - // return true (to zero out static fields - // the very first time // think how to skip recomputing static field in qm scf iterations // check on efp->do_gradient breaks gtests... - if ( if_clean_field(efp, counter_crystal) ) - efp_balance_work(efp, zero_static_field, NULL); - if (res = compute_elec_field_crystal(efp)) return res; diff --git a/src/private.h b/src/private.h index aa8ee10..c848c6e 100644 --- a/src/private.h +++ b/src/private.h @@ -64,12 +64,14 @@ struct polarizable_pt { double x, y, z; mat_t tensor; vec_t elec_field; - vec_t elec_field_wf; + vec_t elec_field_wf; // electric field due to wavefunction of the QM region in the current state vec_t ligand_field; // field due to ligand - vec_t indip; - vec_t indip_old; + vec_t indip; // current induced dipole for the current electronic state + vec_t indip_old; // previous induced dipole for the current electronic state vec_t indipconj; vec_t indipconj_old; + vec_t indip_gs; // induced dipole for the ground electronic state + vec_t indipconj_gs; // induced dipole for the ground electronic state }; /* polarizable point on a ligand to store fields due to other fragments */ @@ -225,12 +227,6 @@ struct efp { /* ligand index in fragment list */ size_t ligand_index; - /* callback which computes electric field from electrons */ - efp_electron_density_field_fn get_electron_density_field; - - /* user data for get_electron_density_field */ - void *get_electron_density_field_user_data; - /* user parameters for this EFP computation */ struct efp_opts opts; diff --git a/tests/symm_2pw_printout.in b/tests/symm_2pw_printout.in new file mode 100644 index 0000000..218b4d0 --- /dev/null +++ b/tests/symm_2pw_printout.in @@ -0,0 +1,337 @@ +run_type sp +ref_energy -2.4694268198 +coord points +terms elec pol disp xr +elec_damp screen +pol_damp tt +disp_damp tt +enable_cutoff true +userlib_path ./ +enable_pbc true +enable_pairwise true +ligand 0 +swf_cutoff 11 +periodic_box 24.0048 27.0356 25.7166 90.0 109.424 90.0 +symmetry true +symm_frag frag + + + fragment XXII_exp_N2_1_L +-6.509027 -0.703872 4.315925 +-5.527225 0.739221 4.530684 +-6.037509 2.064371 4.56864 + + fragment XXII_exp_N2_1_L +-1.630227 -4.083322 1.7473 +-2.612029 -2.640229 1.532541 +-2.101745 -1.315079 1.494585 + + fragment XXII_exp_N2_2_L +6.509027 0.703872 -4.315925 +5.527225 -0.739221 -4.530684 +6.037509 -2.064371 -4.56864 + + fragment XXII_exp_N2_2_L +1.630227 4.083322 -1.7473 +2.612029 2.640229 -1.532541 +2.101745 1.315079 -1.494585 + + fragment XXII_exp_N2_1_L +-10.785134 -0.703872 16.442375 +-9.803332 0.739221 16.657134 +-10.313616 2.064371 16.69509 + + fragment XXII_exp_N2_1_L +-5.906334 -4.083322 13.87375 +-6.888136 -2.640229 13.658991 +-6.377852 -1.315079 13.621035 + + fragment XXII_exp_N2_2_L +2.232919 0.703872 7.810525 +1.251118 -0.739221 7.595766 +1.761401 -2.064371 7.55781 + + fragment XXII_exp_N2_2_L +-2.64588 4.083322 10.37915 +-1.664079 2.640229 10.593909 +-2.174362 1.315079 10.631865 + + fragment XXII_exp_N2_1_L +-6.509027 6.055028 4.315925 +-5.527225 7.498121 4.530684 +-6.037509 8.823271 4.56864 + + fragment XXII_exp_N2_1_L +-1.630227 2.675578 1.7473 +-2.612029 4.118671 1.532541 +-2.101745 5.443821 1.494585 + + fragment XXII_exp_N2_2_L +6.509027 7.462772 -4.315925 +5.527225 6.019679 -4.530684 +6.037509 4.694529 -4.56864 + + fragment XXII_exp_N2_2_L +1.630227 10.842222 -1.7473 +2.612029 9.399129 -1.532541 +2.101745 8.073979 -1.494585 + + fragment XXII_exp_N2_1_L +-10.785134 6.055028 16.442375 +-9.803332 7.498121 16.657134 +-10.313616 8.823271 16.69509 + + fragment XXII_exp_N2_1_L +-5.906334 2.675578 13.87375 +-6.888136 4.118671 13.658991 +-6.377852 5.443821 13.621035 + + fragment XXII_exp_N2_2_L +2.232919 7.462772 7.810525 +1.251118 6.019679 7.595766 +1.761401 4.694529 7.55781 + + fragment XXII_exp_N2_2_L +-2.64588 10.842222 10.37915 +-1.664079 9.399129 10.593909 +-2.174362 8.073979 10.631865 + + fragment XXII_exp_N2_1_L +-6.509027 12.813928 4.315925 +-5.527225 14.257021 4.530684 +-6.037509 15.582171 4.56864 + + fragment XXII_exp_N2_1_L +-1.630227 9.434478 1.7473 +-2.612029 10.877571 1.532541 +-2.101745 12.202721 1.494585 + + fragment XXII_exp_N2_2_L +6.509027 14.221672 -4.315925 +5.527225 12.778579 -4.530684 +6.037509 11.453429 -4.56864 + + fragment XXII_exp_N2_2_L +1.630227 17.601122 -1.7473 +2.612029 16.158029 -1.532541 +2.101745 14.832879 -1.494585 + + fragment XXII_exp_N2_1_L +-10.785134 12.813928 16.442375 +-9.803332 14.257021 16.657134 +-10.313616 15.582171 16.69509 + + fragment XXII_exp_N2_1_L +-5.906334 9.434478 13.87375 +-6.888136 10.877571 13.658991 +-6.377852 12.202721 13.621035 + + fragment XXII_exp_N2_2_L +2.232919 14.221672 7.810525 +1.251118 12.778579 7.595766 +1.761401 11.453429 7.55781 + + fragment XXII_exp_N2_2_L +-2.64588 17.601122 10.37915 +-1.664079 16.158029 10.593909 +-2.174362 14.832879 10.631865 + + fragment XXII_exp_N2_1_L +-6.509027 19.572828 4.315925 +-5.527225 21.015921 4.530684 +-6.037509 22.341071 4.56864 + + fragment XXII_exp_N2_1_L +-1.630227 16.193378 1.7473 +-2.612029 17.636471 1.532541 +-2.101745 18.961621 1.494585 + + fragment XXII_exp_N2_2_L +6.509027 20.980572 -4.315925 +5.527225 19.537479 -4.530684 +6.037509 18.212329 -4.56864 + + fragment XXII_exp_N2_2_L +1.630227 24.360022 -1.7473 +2.612029 22.916929 -1.532541 +2.101745 21.591779 -1.494585 + + fragment XXII_exp_N2_1_L +-10.785134 19.572828 16.442375 +-9.803332 21.015921 16.657134 +-10.313616 22.341071 16.69509 + + fragment XXII_exp_N2_1_L +-5.906334 16.193378 13.87375 +-6.888136 17.636471 13.658991 +-6.377852 18.961621 13.621035 + + fragment XXII_exp_N2_2_L +2.232919 20.980572 7.810525 +1.251118 19.537479 7.595766 +1.761401 18.212329 7.55781 + + fragment XXII_exp_N2_2_L +-2.64588 24.360022 10.37915 +-1.664079 22.916929 10.593909 +-2.174362 21.591779 10.631865 + + fragment XXII_exp_N2_1_L +5.493373 -0.703872 4.315925 +6.475175 0.739221 4.530684 +5.964891 2.064371 4.56864 + + fragment XXII_exp_N2_1_L +10.372173 -4.083322 1.7473 +9.390371 -2.640229 1.532541 +9.900655 -1.315079 1.494585 + + fragment XXII_exp_N2_2_L +18.511427 0.703872 -4.315925 +17.529625 -0.739221 -4.530684 +18.039909 -2.064371 -4.56864 + + fragment XXII_exp_N2_2_L +13.632627 4.083322 -1.7473 +14.614429 2.640229 -1.532541 +14.104145 1.315079 -1.494585 + + fragment XXII_exp_N2_1_L +1.217266 -0.703872 16.442375 +2.199068 0.739221 16.657134 +1.688784 2.064371 16.69509 + + fragment XXII_exp_N2_1_L +6.096066 -4.083322 13.87375 +5.114264 -2.640229 13.658991 +5.624548 -1.315079 13.621035 + + fragment XXII_exp_N2_2_L +14.235319 0.703872 7.810525 +13.253518 -0.739221 7.595766 +13.763801 -2.064371 7.55781 + + fragment XXII_exp_N2_2_L +9.35652 4.083322 10.37915 +10.338321 2.640229 10.593909 +9.828038 1.315079 10.631865 + + fragment XXII_exp_N2_1_L +5.493373 6.055028 4.315925 +6.475175 7.498121 4.530684 +5.964891 8.823271 4.56864 + + fragment XXII_exp_N2_1_L +10.372173 2.675578 1.7473 +9.390371 4.118671 1.532541 +9.900655 5.443821 1.494585 + + fragment XXII_exp_N2_2_L +18.511427 7.462772 -4.315925 +17.529625 6.019679 -4.530684 +18.039909 4.694529 -4.56864 + + fragment XXII_exp_N2_2_L +13.632627 10.842222 -1.7473 +14.614429 9.399129 -1.532541 +14.104145 8.073979 -1.494585 + + fragment XXII_exp_N2_1_L +1.217266 6.055028 16.442375 +2.199068 7.498121 16.657134 +1.688784 8.823271 16.69509 + + fragment XXII_exp_N2_1_L +6.096066 2.675578 13.87375 +5.114264 4.118671 13.658991 +5.624548 5.443821 13.621035 + + fragment XXII_exp_N2_2_L +14.235319 7.462772 7.810525 +13.253518 6.019679 7.595766 +13.763801 4.694529 7.55781 + + fragment XXII_exp_N2_2_L +9.35652 10.842222 10.37915 +10.338321 9.399129 10.593909 +9.828038 8.073979 10.631865 + + fragment XXII_exp_N2_1_L +5.493373 12.813928 4.315925 +6.475175 14.257021 4.530684 +5.964891 15.582171 4.56864 + + fragment XXII_exp_N2_1_L +10.372173 9.434478 1.7473 +9.390371 10.877571 1.532541 +9.900655 12.202721 1.494585 + + fragment XXII_exp_N2_2_L +18.511427 14.221672 -4.315925 +17.529625 12.778579 -4.530684 +18.039909 11.453429 -4.56864 + + fragment XXII_exp_N2_2_L +13.632627 17.601122 -1.7473 +14.614429 16.158029 -1.532541 +14.104145 14.832879 -1.494585 + + fragment XXII_exp_N2_1_L +1.217266 12.813928 16.442375 +2.199068 14.257021 16.657134 +1.688784 15.582171 16.69509 + + fragment XXII_exp_N2_1_L +6.096066 9.434478 13.87375 +5.114264 10.877571 13.658991 +5.624548 12.202721 13.621035 + + fragment XXII_exp_N2_2_L +14.235319 14.221672 7.810525 +13.253518 12.778579 7.595766 +13.763801 11.453429 7.55781 + + fragment XXII_exp_N2_2_L +9.35652 17.601122 10.37915 +10.338321 16.158029 10.593909 +9.828038 14.832879 10.631865 + + fragment XXII_exp_N2_1_L +5.493373 19.572828 4.315925 +6.475175 21.015921 4.530684 +5.964891 22.341071 4.56864 + + fragment XXII_exp_N2_1_L +10.372173 16.193378 1.7473 +9.390371 17.636471 1.532541 +9.900655 18.961621 1.494585 + + fragment XXII_exp_N2_2_L +18.511427 20.980572 -4.315925 +17.529625 19.537479 -4.530684 +18.039909 18.212329 -4.56864 + + fragment XXII_exp_N2_2_L +13.632627 24.360022 -1.7473 +14.614429 22.916929 -1.532541 +14.104145 21.591779 -1.494585 + + fragment XXII_exp_N2_1_L +1.217266 19.572828 16.442375 +2.199068 21.015921 16.657134 +1.688784 22.341071 16.69509 + + fragment XXII_exp_N2_1_L +6.096066 16.193378 13.87375 +5.114264 17.636471 13.658991 +5.624548 18.961621 13.621035 + + fragment XXII_exp_N2_2_L +14.235319 20.980572 7.810525 +13.253518 19.537479 7.595766 +13.763801 18.212329 7.55781 + + fragment XXII_exp_N2_2_L +9.35652 24.360022 10.37915 +10.338321 22.916929 10.593909 +9.828038 21.591779 10.631865