From 263dbdf5b807bb420e7b3139f2adf5e76ee7d803 Mon Sep 17 00:00:00 2001 From: slipchenko Date: Mon, 8 Aug 2022 23:28:02 -0400 Subject: [PATCH 1/2] matched with qchem at Aug. 8 2022 --- efpmd/src/main.c | 16 +- efpmd/src/parse.c | 4 + src/efp.c | 835 +++++++++++++++++++++++++++--------- src/efp.h | 275 ++++++++++-- src/elec.c | 598 ++++++++++++-------------- src/parse.c | 237 ++++++++++- src/pol.c | 851 +++++++++++++++---------------------- src/poldirect.c | 62 ++- src/private.h | 81 ++-- tests/pairwise_0.in | 2 + tests/symm_2pw_printout.in | 337 +++++++++++++++ 11 files changed, 2191 insertions(+), 1107 deletions(-) create mode 100644 tests/symm_2pw_printout.in diff --git a/efpmd/src/main.c b/efpmd/src/main.c index f3f85c6c..8b9795e5 100644 --- a/efpmd/src/main.c +++ b/efpmd/src/main.c @@ -152,7 +152,7 @@ static struct cfg *make_cfg(void) cfg_add_double(cfg, "thermostat_tau", 1.0e3); cfg_add_double(cfg, "barostat_tau", 1.0e4); - cfg_add_int(cfg, "ligand", 0); + cfg_add_int(cfg, "ligand", -100); cfg_add_bool(cfg, "enable_pairwise", false); cfg_add_bool(cfg, "print_pbc", false); cfg_add_bool(cfg, "symmetry", false); @@ -166,6 +166,8 @@ static struct cfg *make_cfg(void) cfg_add_int(cfg, "update_params", 0); cfg_add_double(cfg, "update_params_cutoff", 0.0); + cfg_add_int(cfg, "print", 0); + return cfg; } @@ -286,12 +288,13 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys) .symmetry = cfg_get_bool(cfg, "symmetry"), .symm_frag = cfg_get_enum(cfg, "symm_frag"), .update_params = cfg_get_int(cfg, "update_params"), - .update_params_cutoff = cfg_get_double(cfg, "update_params_cutoff") + .update_params_cutoff = cfg_get_double(cfg, "update_params_cutoff"), + .print = cfg_get_int(cfg, "print") }; if (opts.xr_cutoff == 0.0) { opts.xr_cutoff = opts.swf_cutoff; - printf("xr_cutoff is set to %lf \n\n\n", opts.xr_cutoff*0.52917721092); + printf("xr_cutoff is set to %lf \n\n", opts.xr_cutoff*0.52917721092); } enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord"); @@ -306,7 +309,7 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys) add_potentials(efp, cfg, sys); for (size_t i = 0; i < sys->n_frags; i++) - check_fail(efp_add_fragment(efp, sys->frags[i].name)); + check_fail(efp_add_fragment(efp, sys->frags[i].name)); if (sys->n_charges > 0) { double q[sys->n_charges]; @@ -331,7 +334,10 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys) if (cfg_get_bool(cfg, "enable_ff")) opts.terms &= ~(EFP_TERM_ELEC | EFP_TERM_POL | EFP_TERM_DISP | EFP_TERM_XR); - check_fail(efp_set_opts(efp, &opts)); + if (opts.enable_pairwise) + check_fail(efp_add_ligand(efp, opts.ligand)); + + check_fail(efp_set_opts(efp, &opts)); check_fail(efp_prepare(efp)); check_fail(efp_set_symmlist(efp)); diff --git a/efpmd/src/parse.c b/efpmd/src/parse.c index 5b81dc1f..f996db11 100644 --- a/efpmd/src/parse.c +++ b/efpmd/src/parse.c @@ -253,6 +253,10 @@ struct sys *parse_input(struct cfg *cfg, const char *path) if (cfg_get_bool(cfg, "enable_pbc")) cfg_set_bool(cfg, "enable_cutoff", true); + // turn off pairwise calculations if ligand is not set + if (cfg_get_int(cfg, "ligand") == -100 && cfg_get_bool(cfg, "enable_pairwise") == true) + error("Specify ligand for pairwise calculations"); + check_cfg(cfg); efp_stream_close(stream); diff --git a/src/efp.c b/src/efp.c index 80ee95f9..a20ac596 100644 --- a/src/efp.c +++ b/src/efp.c @@ -35,20 +35,6 @@ #include "stream.h" #include "util.h" -static void -add_screen2_params(struct frag *frag) { - double *scr; - scr = (double *)malloc(frag->n_multipole_pts * sizeof(double)); - // if (scr == NULL) - // return EFP_RESULT_NO_MEMORY; - for (int i=0; in_multipole_pts; i++){ - scr[i] = 10.0; // assign large value - no effective screening - } - - if (frag->screen_params) - free(frag->screen_params); - frag->screen_params = scr; -} static void update_fragment(struct frag *frag) @@ -194,8 +180,6 @@ free_frag(struct frag *frag) free(frag->xr_fock_mat); free(frag->xr_wf); free(frag->xrfit); - free(frag->screen_params); - free(frag->ai_screen_params); for (size_t i = 0; i < 3; i++) free(frag->xr_wf_deriv[i]); @@ -211,6 +195,21 @@ free_frag(struct frag *frag) /* don't do free(frag) here */ } +static void +free_ligand(struct ligand *ligand) +{ + if (!ligand) + return; + + //free_frag(&ligand->ligand_frag); + for (size_t i=0; in_ligand_pts; i++) { + if (ligand->ligand_pts[i].fragment_field) + free(ligand->ligand_pts[i].fragment_field); + } + if (ligand->ligand_pts) + free(ligand->ligand_pts); +} + static enum efp_result copy_frag(struct frag *dest, const struct frag *src) { @@ -232,20 +231,6 @@ copy_frag(struct frag *dest, const struct frag *src) return EFP_RESULT_NO_MEMORY; memcpy(dest->multipole_pts, src->multipole_pts, size); } - if (src->screen_params) { - size = src->n_multipole_pts * sizeof(double); - dest->screen_params = (double *)malloc(size); - if (!dest->screen_params) - return EFP_RESULT_NO_MEMORY; - memcpy(dest->screen_params, src->screen_params, size); - } - if (src->ai_screen_params) { - size = src->n_multipole_pts * sizeof(double); - dest->ai_screen_params = (double *)malloc(size); - if (!dest->ai_screen_params) - return EFP_RESULT_NO_MEMORY; - memcpy(dest->ai_screen_params, src->ai_screen_params, size); - } if (src->polarizable_pts) { size = src->n_polarizable_pts * sizeof(struct polarizable_pt); dest->polarizable_pts = (struct polarizable_pt *)malloc(size); @@ -333,6 +318,36 @@ copy_frag(struct frag *dest, const struct frag *src) return EFP_RESULT_SUCCESS; } +/* +static enum efp_result +copy_ligand(struct ligand *dest, const struct ligand *src) { + size_t size; + + memcpy(dest, src, sizeof(*dest)); + + if (src->ligand_frag) + copy_frag(dest->ligand_frag, src->ligand_frag); + + if (src->ligand_pts) { + size = src->n_ligand_pts * sizeof(struct ligand_pt); + dest->ligand_pts = (struct ligand_pt *) malloc(size); + if (!dest->ligand_pts) + return EFP_RESULT_NO_MEMORY; + memcpy(dest->ligand_pts, src->ligand_pts, size); + + for (size_t i=0; in_ligand_pts; i++) { + const struct ligand_pt *src_pt = src->ligand_pts + i; + struct ligand_pt *dest_pt = dest->ligand_pts + i; + size = src_pt->n_frag * sizeof(vec_t); + dest_pt->fragment_field = (vec_t *)malloc(size); + if (!dest_pt->fragment_field) + return EFP_RESULT_NO_MEMORY; + memcpy(dest_pt->fragment_field, src_pt->fragment_field, size); + } + } +} +*/ + // updates (shifts) parameters of fragment based on coordinates of fragment atoms static enum efp_result update_params(struct efp_atom *atoms, const struct frag *lib_orig, const struct frag *lib_current) { @@ -394,31 +409,35 @@ check_frag_params(const struct efp_opts *opts, struct frag *frag) { if ((opts->terms & EFP_TERM_ELEC) || (opts->terms & EFP_TERM_AI_ELEC)) { if (!frag->multipole_pts) { - efp_log("electrostatic parameters are missing"); - return EFP_RESULT_FATAL; - } - if (opts->elec_damp == EFP_ELEC_DAMP_SCREEN && - frag->screen_params == NULL) { - efp_log("electrostatic screening parameters are missing; continue"); - add_screen2_params(frag); - //return EFP_RESULT_FATAL; + printf("WARNING! Multipole parameters are missing on fragment %s\n", frag->name); + // efp_log("electrostatic parameters are missing"); + // return EFP_RESULT_FATAL; } } if ((opts->terms & EFP_TERM_POL) || (opts->terms & EFP_TERM_AI_POL)) { - if (!frag->polarizable_pts || !frag->multipole_pts) { - efp_log("polarization parameters are missing"); - return EFP_RESULT_FATAL; + if (!frag->multipole_pts) + printf("WARNING! Multipole parameters are missing on fragment %s\n", frag->name); + if (!frag->polarizable_pts) { + printf("WARNING! Polarizability parameters are missing on fragment %s\n", frag->name); + + // efp_log("polarization parameters are missing"); + // return EFP_RESULT_FATAL; } } if ((opts->terms & EFP_TERM_DISP) || (opts->terms & EFP_TERM_AI_DISP)) { if (frag->dynamic_polarizable_pts == NULL) { - efp_log("dispersion parameters are missing"); - return EFP_RESULT_FATAL; + printf("WARNING! Dynamic polarizability parameters are " + "missing on fragment %s\n", frag->name); + + // efp_log("dispersion parameters are missing"); + // return EFP_RESULT_FATAL; } if (opts->disp_damp == EFP_DISP_DAMP_OVERLAP && frag->n_lmo != frag->n_dynamic_polarizable_pts) { - efp_log("number of polarization points does not " - "match number of LMOs"); + //efp_log("number of polarization points does not " + // "match number of LMOs"); + printf("FATAL ERROR! The number of polarization points does not " + "match the number of LMOs of fragment %s\n", frag->name); return EFP_RESULT_FATAL; } } @@ -427,8 +446,10 @@ check_frag_params(const struct efp_opts *opts, struct frag *frag) !frag->xr_fock_mat || !frag->xr_wf || !frag->lmo_centroids) { - efp_log("exchange repulsion parameters are missing"); - return EFP_RESULT_FATAL; + printf("WARNING! Exchange-repulsion parameters are " + "missing on fragment %s\n", frag->name); + //efp_log("exchange repulsion parameters are missing"); + //return EFP_RESULT_FATAL; } } return EFP_RESULT_SUCCESS; @@ -439,12 +460,14 @@ check_params(struct efp *efp) { enum efp_result res; - for (size_t i = 0; i < efp->n_frag; i++) - if ((res = check_frag_params(&efp->opts, efp->frags + i))) { + for (size_t i = 0; i < efp->n_frag; i++) { + if (efp->opts.print > 1) + print_frag_info(efp, i); + if ((res = check_frag_params(&efp->opts, efp->frags + i))) { efp_log("check_params() failure"); return res; - } - + } + } return EFP_RESULT_SUCCESS; } @@ -479,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 @@ -496,55 +522,61 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, size_t n_lmo_ij = efp->frags[i].n_lmo * efp->frags[fr_j].n_lmo; - s = (double *)calloc(n_lmo_ij, sizeof(double)); - ds = (six_t *)calloc(n_lmo_ij, sizeof(six_t)); - - if (do_xr(&efp->opts)) { - double exr, ecp; - - efp_frag_frag_xr(efp, i, fr_j, - s, ds, &exr, &ecp); - e_xr += exr; - e_cp += ecp; - - /* */ - if (efp->opts.enable_pairwise) { - if (i == efp->opts.ligand) { - efp->pair_energies[fr_j].exchange_repulsion = exr; - efp->pair_energies[fr_j].charge_penetration = ecp; - } - if (fr_j == efp->opts.ligand) { - efp->pair_energies[i].exchange_repulsion = exr; - efp->pair_energies[i].charge_penetration = ecp; + // ugly but simple check when not to compute exchange-repulsion and overlap screening + if (n_lmo_ij > 0) { + s = (double *) calloc(n_lmo_ij, sizeof(double)); + ds = (six_t *) calloc(n_lmo_ij, sizeof(six_t)); + + if (do_xr(&efp->opts)) { + double exr, ecp; + + efp_frag_frag_xr(efp, i, fr_j, + s, ds, &exr, &ecp); + e_xr += exr; + e_cp += ecp; + + /* */ + if (if_pairwise) { + if (i == efp->ligand_index) { + efp->pair_energies[fr_j].exchange_repulsion = exr; + efp->pair_energies[fr_j].charge_penetration = ecp; + } + if (fr_j == efp->ligand_index) { + efp->pair_energies[i].exchange_repulsion = exr; + efp->pair_energies[i].charge_penetration = ecp; + } } } - } - if (do_elec(&efp->opts)) { - e_elec_tmp = efp_frag_frag_elec(efp, - i, fr_j); + } + 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 += e_elec_tmp; /* */ - if (efp->opts.enable_pairwise) { - if (i == efp->opts.ligand) + if (if_pairwise) { + if (i == efp->ligand_index) efp->pair_energies[fr_j].electrostatic = e_elec_tmp; - if (fr_j == efp->opts.ligand) + if (fr_j == efp->ligand_index) efp->pair_energies[i].electrostatic = e_elec_tmp; } } - if (do_disp(&efp->opts)) { + if (do_disp(&efp->opts) && efp->frags[i].n_dynamic_polarizable_pts > 0 && + efp->frags[fr_j].n_dynamic_polarizable_pts > 0) { e_disp_tmp = efp_frag_frag_disp(efp, 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) efp->pair_energies[i].dispersion = e_disp_tmp; } } - free(s); - free(ds); + if (n_lmo_ij > 0) { + free(s); + free(ds); + } } } } @@ -552,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 @@ -582,33 +620,35 @@ compute_two_body_crystal(struct efp *efp) six_t *ds; size_t n_lmo_ij = efp->frags[i].n_lmo * efp->frags[fr_j].n_lmo; - - s = (double *)calloc(n_lmo_ij, sizeof(double)); - ds = (six_t *)calloc(n_lmo_ij, sizeof(six_t)); - - if (do_xr(&efp->opts)) { - double exr, ecp; - - efp_frag_frag_xr(efp, i, fr_j, - s, ds, &exr, &ecp); - e_xr += exr * factor; - e_cp += ecp * factor; - - /* */ - if (efp->opts.enable_pairwise) { - if (i == efp->opts.ligand) { - efp->pair_energies[fr_j].exchange_repulsion = exr; - efp->pair_energies[fr_j].charge_penetration = ecp; - } - if (fr_j == efp->opts.ligand) { - efp->pair_energies[i].exchange_repulsion = exr; - efp->pair_energies[i].charge_penetration = ecp; + // if need to compute exrep and overlap screening for this pair + if (n_lmo_ij > 0) { + s = (double *) calloc(n_lmo_ij, sizeof(double)); + ds = (six_t *) calloc(n_lmo_ij, sizeof(six_t)); + + if (do_xr(&efp->opts)) { + double exr, ecp; + + efp_frag_frag_xr(efp, i, fr_j, + s, ds, &exr, &ecp); + e_xr += exr * factor; + e_cp += ecp * factor; + + /* */ + if (efp->opts.enable_pairwise) { + if (i == efp->opts.ligand) { + efp->pair_energies[fr_j].exchange_repulsion = exr; + efp->pair_energies[fr_j].charge_penetration = ecp; + } + if (fr_j == efp->opts.ligand) { + efp->pair_energies[i].exchange_repulsion = exr; + efp->pair_energies[i].charge_penetration = ecp; + } } } } - if (do_elec(&efp->opts)) { - e_elec_tmp = efp_frag_frag_elec(efp, - i, fr_j); + 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 += e_elec_tmp * factor; /* */ @@ -619,9 +659,9 @@ compute_two_body_crystal(struct efp *efp) efp->pair_energies[i].electrostatic = e_elec_tmp; } } - if (do_disp(&efp->opts)) { - e_disp_tmp = efp_frag_frag_disp(efp, - i, fr_j, s, ds); + if (do_disp(&efp->opts) && efp->frags[i].n_dynamic_polarizable_pts > 0 && + efp->frags[fr_j].n_dynamic_polarizable_pts > 0) { + e_disp_tmp = efp_frag_frag_disp(efp, i, fr_j, s, ds); e_disp += e_disp_tmp * factor; /* */ if (efp->opts.enable_pairwise) { @@ -631,9 +671,10 @@ compute_two_body_crystal(struct efp *efp) efp->pair_energies[i].dispersion = e_disp_tmp; } } - free(s); - free(ds); - + if (n_lmo_ij > 0) { + free(s); + free(ds); + } } } } @@ -1101,10 +1142,9 @@ efp_get_stress_tensor(struct efp *efp, double *stress) } EFP_EXPORT enum efp_result -efp_get_ai_screen(struct efp *efp, size_t frag_idx, double *screen) +efp_get_frag_ai_screen(struct efp *efp, size_t frag_idx, double *screen, int if_screen) { const struct frag *frag; - size_t size; assert(efp); assert(screen); @@ -1112,43 +1152,70 @@ efp_get_ai_screen(struct efp *efp, size_t frag_idx, double *screen) frag = &efp->frags[frag_idx]; - if (frag->ai_screen_params == NULL) { - efp_log("no screening parameters found for %s", frag->name); - return EFP_RESULT_FATAL; - } + if_screen = 0; + for (int i=0; in_multipole_pts; i++) { + struct multipole_pt *pt = frag->multipole_pts + i; - size = frag->n_multipole_pts * sizeof(double); - memcpy(screen, frag->ai_screen_params, size); + *screen++ = pt->screen0; + if (pt->if_scr0) + if_screen = 1; + } return EFP_RESULT_SUCCESS; } +EFP_EXPORT enum efp_result +efp_get_all_ai_screen(struct efp *efp, double *screen) +{ + assert(efp); + assert(screen); + + 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; + *screen++ = pt->screen0; + } + } + + return EFP_RESULT_SUCCESS; +} +/* +EFP_EXPORT enum efp_result +efp_get_ai_screen(struct efp *efp, size_t frag_idx, double *screen) +{ + const struct frag *frag; + size_t size; + + assert(efp); + assert(screen); + assert(frag_idx < efp->n_frag); + + frag = &efp->frags[frag_idx]; + + if (frag->ai_screen_params == NULL) { + efp_log("no screening parameters found for %s", frag->name); + return EFP_RESULT_FATAL; + } + + size = frag->n_multipole_pts * sizeof(double); + memcpy(screen, frag->ai_screen_params, size); + + return EFP_RESULT_SUCCESS; +} +*/ EFP_EXPORT enum efp_result efp_prepare(struct efp *efp) { assert(efp); efp->n_polarizable_pts = 0; - for (size_t i = 0; i < efp->n_frag; i++) { efp->frags[i].polarizable_offset = efp->n_polarizable_pts; 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->indip = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); - efp->indipconj = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); - efp->indip_old = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); - efp->indipconj_old = (vec_t *)calloc(efp->n_polarizable_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)); @@ -1235,6 +1302,8 @@ efp_compute(struct efp *efp, int do_gradient) assert(efp); + static int efp_counter = 0; + if (efp->grad == NULL) { efp_log("call efp_prepare after all fragments are added"); return EFP_RESULT_FATAL; @@ -1242,10 +1311,11 @@ efp_compute(struct efp *efp, int do_gradient) efp->do_gradient = do_gradient; - if ((res = check_params(efp))) { - efp_log("check_params() failure"); - return res; - } + if (efp_counter == 0) + if ((res = check_params(efp))) { + efp_log("check_params() failure"); + return res; + } memset(&efp->energy, 0, sizeof(efp->energy)); memset(&efp->stress, 0, sizeof(efp->stress)); @@ -1296,6 +1366,8 @@ efp_compute(struct efp *efp, int do_gradient) efp->energy.ai_dispersion + efp->energy.exchange_repulsion; + efp_counter++; + return EFP_RESULT_SUCCESS; } @@ -1344,38 +1416,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; @@ -1441,6 +1503,62 @@ efp_get_multipole_values(struct efp *efp, double *mult) return EFP_RESULT_SUCCESS; } +EFP_EXPORT enum efp_result +efp_get_dipole_values(struct efp *efp, double *dipoles) +{ + assert(efp); + assert(dipoles); + + 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; + + *dipoles++ = pt->dipole.x; + *dipoles++ = pt->dipole.y; + *dipoles++ = pt->dipole.z; + } + } + return EFP_RESULT_SUCCESS; +} + +EFP_EXPORT enum efp_result +efp_get_quadrupole_values(struct efp *efp, double *quad) +{ + assert(efp); + assert(quad); + + 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; + for (size_t t = 0; t < 6; t++) + *quad++ = pt->quadrupole[t]; + } + } + return EFP_RESULT_SUCCESS; +} + +EFP_EXPORT enum efp_result +efp_get_octupole_values(struct efp *efp, double *oct) +{ + assert(efp); + assert(oct); + + 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; + for (size_t t = 0; t < 10; t++) + *oct++ = pt->octupole[t]; + } + } + return EFP_RESULT_SUCCESS; +} + EFP_EXPORT enum efp_result efp_get_frag_induced_dipole_count(struct efp *efp, size_t frag_idx, size_t *n_dip) { @@ -1491,41 +1609,109 @@ efp_get_induced_dipole_coordinates(struct efp *efp, double *xyz) EFP_EXPORT enum efp_result efp_get_induced_dipole_values(struct efp *efp, double *dip) { - assert(efp); - assert(dip); + assert(efp); + assert(dip); - memcpy(dip, efp->indip, efp->n_polarizable_pts * sizeof(vec_t)); + 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; + + *dip++ = pt->indip.x; + *dip++ = pt->indip.y; + *dip++ = pt->indip.z; + } + } return EFP_RESULT_SUCCESS; } EFP_EXPORT enum efp_result efp_get_induced_dipole_conj_values(struct efp *efp, double *dip) { - assert(efp); - assert(dip); + assert(efp); + assert(dip); + + for (size_t i = 0; i < efp->n_frag; i++) { + struct frag *frag = efp->frags + i; - memcpy(dip, efp->indipconj, efp->n_polarizable_pts * sizeof(vec_t)); + for (size_t j = 0; j < frag->n_polarizable_pts; j++) { + struct polarizable_pt *pt = frag->polarizable_pts + j; + + *dip++ = pt->indipconj.x; + *dip++ = pt->indipconj.y; + *dip++ = pt->indipconj.z; + } + } return EFP_RESULT_SUCCESS; } EFP_EXPORT enum efp_result -efp_get_old_induced_dipole_values(struct efp *efp, double *dip) +efp_set_induced_dipole_values(struct efp *efp, double *dip, int if_conjug) { assert(efp); assert(dip); - memcpy(dip, efp->indip_old, efp->n_polarizable_pts * sizeof(vec_t)); + size_t index = 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; + + if (if_conjug) { + pt->indipconj.x = dip[index]; + pt->indipconj.y = dip[index + 1]; + pt->indipconj.z = dip[index + 2]; + } + else { + pt->indip.x = dip[index]; + pt->indip.y = dip[index + 1]; + pt->indip.z = dip[index + 2]; + } + index+=3; + } + } return EFP_RESULT_SUCCESS; } + +EFP_EXPORT enum efp_result +efp_get_old_induced_dipole_values(struct efp *efp, double *dip) +{ + assert(efp); + assert(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; + + *dip++ = pt->indip_old.x; + *dip++ = pt->indip_old.y; + *dip++ = pt->indip_old.z; + } + } +} + EFP_EXPORT enum efp_result efp_get_old_induced_dipole_conj_values(struct efp *efp, double *dip) { assert(efp); assert(dip); - memcpy(dip, efp->indipconj_old, efp->n_polarizable_pts * sizeof(vec_t)); - return EFP_RESULT_SUCCESS; + 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; + + *dip++ = pt->indipconj_old.x; + *dip++ = pt->indipconj_old.y; + *dip++ = pt->indipconj_old.z; + } + } } EFP_EXPORT enum efp_result @@ -1593,6 +1779,9 @@ efp_shutdown(struct efp *efp) free_frag(efp->lib_current[i]); free(efp->lib_current[i]); } + + free_ligand(efp->ligand); + free(efp->ligand); free(efp->frags); free(efp->lib); free(efp->lib_current); @@ -1600,14 +1789,9 @@ efp_shutdown(struct efp *efp) free(efp->ptc); free(efp->ptc_xyz); free(efp->ptc_grad); - free(efp->indip); - free(efp->indipconj); - free(efp->indip_old); - free(efp->indipconj_old); 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); @@ -1710,6 +1894,47 @@ efp_add_fragment(struct efp *efp, const char *name) return EFP_RESULT_SUCCESS; } + +EFP_EXPORT enum efp_result +efp_add_ligand(struct efp *efp, int ligand_index) { + + assert(efp); + assert(ligand_index >= -1); + + // ligand_index=-1 means ligand is QM; -100 is default when ligand is not specified + // 0 and larger: ligand is a fragment with with index + if (ligand_index > -1) { + assert( (size_t)ligand_index < efp->n_frag); + efp->ligand_index = (size_t)ligand_index; + + efp->ligand = (struct ligand *) calloc(1, sizeof(struct ligand)); + if (efp->ligand == NULL) + return EFP_RESULT_NO_MEMORY; + + struct ligand *lig = efp->ligand; + + // copy_frag(lig->ligand_frag, &efp->frags[ligand_index]); + lig->ligand_frag = &efp->frags[ligand_index]; + lig->n_ligand_pts = efp->frags[ligand_index].n_polarizable_pts; + + size_t size; + size = lig->n_ligand_pts * sizeof(struct ligand_pt); + lig->ligand_pts = (struct ligand_pt *) malloc(size); + if (lig->ligand_pts == NULL) + return EFP_RESULT_NO_MEMORY; + + for (size_t i = 0; i < lig->n_ligand_pts; i++) { + struct ligand_pt *pt = lig->ligand_pts + i; + size = efp->n_frag * sizeof(vec_t); + pt->n_frag = efp->n_frag; + pt->fragment_field = (vec_t *) malloc(size); + if (!pt->fragment_field) + return EFP_RESULT_NO_MEMORY; + } + } + return EFP_RESULT_SUCCESS; +} + EFP_EXPORT enum efp_result efp_skip_fragments(struct efp *efp, size_t i, size_t j, int value) { @@ -1737,27 +1962,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) { @@ -1873,6 +2077,104 @@ efp_set_frag_atoms(struct efp *efp, size_t frag_idx, size_t n_atoms, return EFP_RESULT_SUCCESS; } +EFP_EXPORT 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) +{ + assert(efp); + assert(mult_pt); + assert(frag_idx < efp->n_frag); + assert(pt_idx < efp->frags[frag_idx].n_multipole_pts); + + struct frag *frag; + frag = efp->frags + frag_idx; + struct multipole_pt *pt; + pt = frag->multipole_pts + pt_idx; + + mult_pt->x = pt->x; + mult_pt->y = pt->y; + mult_pt->z = pt->z; + mult_pt->znuc = pt->znuc; + mult_pt->monopole = pt->monopole; + mult_pt->dipole[0] = pt->dipole.x; + mult_pt->dipole[1] = pt->dipole.y; + mult_pt->dipole[2] = pt->dipole.z; + for (size_t i=0; i<6; i++) + mult_pt->quadrupole[i] = pt->quadrupole[i]; + for (size_t i=0; i<10; i++) + mult_pt->octupole[i] = pt->octupole[i]; + mult_pt->screen0 = pt->screen0; + mult_pt->if_screen = pt->if_scr0; + mult_pt->rank = 0; + if ( pt->if_dip ) + mult_pt->rank = 1; + if ( pt->if_quad ) + mult_pt->rank = 2; + if ( pt->if_oct ) + mult_pt->rank = 3; + + 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) @@ -1901,6 +2203,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" "\n" "Journal References:\n" " - Kaliman and Slipchenko, JCC 2013.\n" @@ -2090,4 +2393,136 @@ n_symm_frag(struct efp *efp, size_t *symm_frag) { } } +void +print_mult_pt(struct efp *efp, size_t frag_index, size_t pt_index) { + struct multipole_pt *pt = efp->frags[frag_index].multipole_pts + pt_index; + printf(" Multipole point %s %lf %lf %lf\n", pt->label, pt->x, pt->y, pt->z); + printf(" znuc, monopole %lf %lf\n", pt->znuc, pt->monopole); + printf(" dipole %lf %lf %lf\n", pt->dipole.x, pt->dipole.y, pt->dipole.z); + printf(" duadrupole %lf %lf %lf %lf %lf %lf\n", pt->quadrupole[0], pt->quadrupole[1], + pt->quadrupole[2], pt->quadrupole[3], pt->quadrupole[4], pt->quadrupole[5]); + printf(" octupole %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + pt->octupole[0], pt->octupole[1], pt->octupole[2], pt->octupole[3], pt->octupole[4], + pt->octupole[5], pt->octupole[6], pt->octupole[7], pt->octupole[8], pt->octupole[9]); + printf(" screen, screen2 %lf %lf\n", pt->screen0, pt->screen2); + printf("\n"); +} + +void +print_atoms(struct efp *efp, size_t frag_index, size_t atom_index) { + struct efp_atom *atom = efp->frags[frag_index].atoms + atom_index; + printf(" Atom %s %lf %lf %lf\n", atom->label, atom->x, atom->y, atom->z); + printf(" znuc, mass %lf %lf\n", atom->znuc, atom->mass); + printf("\n"); +} + +void +print_pol_pt(struct efp *efp, size_t frag_index, size_t pol_index) { + struct polarizable_pt *pt = efp->frags[frag_index].polarizable_pts + pol_index; + printf(" Polarizability point coords %lf %lf %lf\n", pt->x, pt->y, pt->z); + printf(" polarizability tensor %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + pt->tensor.xx, pt->tensor.xy, pt->tensor.xz, pt->tensor.yx, pt->tensor.yy, pt->tensor.yz, + pt->tensor.zx, pt->tensor.zy, pt->tensor.zz); + printf(" electric field %lf %lf %lf\n", pt->elec_field.x, pt->elec_field.y, pt->elec_field.z); + printf(" wf electric field %lf %lf %lf\n", pt->elec_field_wf.x, pt->elec_field_wf.y, pt->elec_field_wf.z); + printf(" ligand electric field %lf %lf %lf\n", pt->ligand_field.x, pt->ligand_field.y, pt->ligand_field.z); + printf(" induced dipole %lf %lf %lf\n", pt->indip.x, pt->indip.y, pt->indip.z); + 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"); +} + +void print_ligand(struct efp *efp, size_t frag_index) { + if (! efp->opts.enable_pairwise) + return; + if (efp->ligand_index == frag_index) { + printf("Ligand index %d\n", efp->ligand_index); + if (efp->ligand_index > -1) { + printf("Number of ligand points %d", efp->ligand->n_ligand_pts); + } + } +} +void +print_frag_info(struct efp *efp, size_t frag_index) { + struct frag *fr = efp->frags + frag_index; + printf("Fragment %s\n", fr->name); + + for (int i=0; i < fr->n_atoms; i++) { + print_atoms(efp, frag_index, i); + } + + for (int i=0; i < fr->n_multipole_pts; i++) { + print_mult_pt(efp, frag_index, i); + } + + for (int i=0; i < fr->n_polarizable_pts; i++) { + print_pol_pt(efp, frag_index, i); + } + + print_ligand(efp, frag_index); + + printf("\n"); +} + +void +print_efp_mult_pt(struct efp_mult_pt *pt) { + printf(" Multipole point of rank %d\n", pt->rank); + printf(" Coordinates %lf %lf %lf\n", pt->x, pt->y, pt->z); + printf(" znuc, monopole %lf %lf\n", pt->znuc, pt->monopole); + printf(" dipole %lf %lf %lf\n", pt->dipole[0], pt->dipole[1], pt->dipole[2]); + printf(" duadrupole %lf %lf %lf %lf %lf %lf\n", pt->quadrupole[0], pt->quadrupole[1], + pt->quadrupole[2], pt->quadrupole[3], pt->quadrupole[4], pt->quadrupole[5]); + printf(" octupole %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + pt->octupole[0], pt->octupole[1], pt->octupole[2], pt->octupole[3], pt->octupole[4], + pt->octupole[5], pt->octupole[6], pt->octupole[7], pt->octupole[8], pt->octupole[9]); + 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 373147bf..5c55c156 100644 --- a/src/efp.h +++ b/src/efp.h @@ -28,6 +28,7 @@ #define LIBEFP_EFP_H #include +#include /** \file efp.h * Public libefp interface. @@ -177,6 +178,8 @@ struct efp_opts { int update_params; /** Cutoff when updating parameters is "safe". Default 0.0 (never safe) */ double update_params_cutoff; + /** Level of print out */ + int print; }; /** EFP energy terms. */ @@ -224,12 +227,39 @@ struct efp_energy { /** EFP atom info. */ struct efp_atom { - char label[32]; /**< Atom label. */ - double x; /**< X coordinate of atom position. */ - double y; /**< Y coordinate of atom position. */ - double z; /**< Z coordinate of atom position. */ - double mass; /**< Atom mass. */ - double znuc; /**< Nuclear charge. */ + char label[32]; /**< Atom label. */ + double x; /**< X coordinate of atom position. */ + double y; /**< Y coordinate of atom position. */ + double z; /**< Z coordinate of atom position. */ + double mass; /**< Atom mass. */ + double znuc; /**< Nuclear charge. */ +}; + +/** Multipole point for working with external programs */ +struct efp_mult_pt { + double x; /**< X coordinate */ + double y; /**< Y coordinate */ + double z; /**< Z coordinate */ + double znuc; /**< Nuclear charge */ + double monopole; /**< Monopole */ + double dipole[3]; /**< Dipole */ + double quadrupole[6]; /**< Quadrupole */ + double octupole[10]; /**< Octupole */ + size_t rank; /** < Highest non-zero multipole: 0 - monopole, 1 - dipole, 2 - quad, 3 - oct */ + double screen0; /**< AI-EFP screening parameter */ + 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 */ }; /** @@ -254,8 +284,8 @@ struct efp_atom { * \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. @@ -345,6 +375,14 @@ enum efp_result efp_add_potential(struct efp *efp, const char *path); */ enum efp_result efp_add_fragment(struct efp *efp, const char *name); +/** + * Add a ligand fragment to teh system + * @param[in] efp + * @param[in] ligand_index Index of the ligand in the fragment list + * @return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result efp_add_ligand(struct efp *efp, int ligand_index); + /** * Prepare the calculation. * @@ -382,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. @@ -395,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. @@ -646,8 +684,40 @@ enum efp_result efp_get_stress_tensor(struct efp *efp, double *stress); * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -enum efp_result efp_get_ai_screen(struct efp *efp, size_t frag_idx, - double *screen); +//enum efp_result efp_get_ai_screen(struct efp *efp, size_t frag_idx, +// double *screen); + + +/** + * Get all ab initio screening parameters. + * + * \param[in] efp The efp structure. + * + * \param[out] screen Array of N elements where screening parameters will be + * stored. N is the total number of multipole points in all fragments. + * + * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result efp_get_all_ai_screen(struct efp *efp, double *screen); + +/** + * Get the ab initio screening parameters for one 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[out] screen Array of N elements where screening parameters will be + * stored. N is the total number of multipole points in fragment frag_idx. + * + * \param[out] if_screen 0 if screening parameters are set to 10.0 (could be ignored); + * 1 if at least one point has != 10 screening parameter + * + * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result efp_get_frag_ai_screen(struct efp *efp, size_t frag_idx, + double *screen, int if_screen); /** * Set ab initio orbital energies. @@ -767,18 +837,19 @@ enum efp_result efp_get_frag_multiplicity(struct efp *efp, size_t frag_idx, enum efp_result efp_get_frag_multipole_count(struct efp *efp, size_t frag_idx, size_t *n_mult); +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); +enum efp_result +efp_get_frag_rank(struct efp *efp, size_t frag_idx, size_t *rank); /** * Get total number of multipoles from EFP electrostatics. @@ -828,6 +899,55 @@ enum efp_result efp_get_multipole_coordinates(struct efp *efp, double *xyz); */ enum efp_result efp_get_multipole_values(struct efp *efp, double *mult); +/** + * 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 + * 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. + */ +enum efp_result efp_get_dipole_values(struct efp *efp, double *dipoles); + +/** + * Get electrostatics quadrupoles from EFP fragments. + * + * \param[in] efp The efp structure. + * + * \param[out] quad Array with all efp quadrupoles. + * + * The size of the \p mult array must be at least [6 * \p n_mult] elements + * where \p n_mult is the value returned by the ::efp_get_multipole_count function. + * + * * Quadrupoles are stored in the following order: + * \a xx, \a yy, \a zz, \a xy, \a xz, \a yz + * + * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result efp_get_quadrupole_values(struct efp *efp, double *quad); + +/** + * Get electrostatics octupoles from EFP fragments. + * + * \param[in] efp The efp structure. + * + * \param[out] oct Array with all efp octupoles. + * + * The size of the \p mult array must be at least [10 * \p n_mult] elements + * where \p n_mult is the value returned by the ::efp_get_multipole_count function. + * + * Octupoles are stored in the following order: + * \a xxx, \a yyy, \a zzz, \a xxy, \a xxz, + * \a xyy, \a yyz, \a xzz, \a yzz, \a xyz + * + * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +enum efp_result efp_get_octupole_values(struct efp *efp, double *oct); + /** * Get the number of polarization induced dipoles from a particular fragment. * @@ -910,6 +1030,15 @@ enum efp_result efp_get_old_induced_dipole_values(struct efp *efp, double *dip); */ enum efp_result efp_get_old_induced_dipole_conj_values(struct efp *efp, double *dip); +/** + * Writes induced dipoles from dip array into polarizable points + * @param efp + * @param dip pointer to array with induced dipoles + * @param if_conjug 0 to write indip and 1 to write indipconj + * @return + */ +enum efp_result efp_set_induced_dipole_values(struct efp *efp, double *dip, int if_conjug); + /** * Get the number of LMOs in a fragment. * @@ -1108,6 +1237,42 @@ enum efp_result efp_get_frag_atoms(struct efp *efp, size_t frag_idx, enum efp_result efp_set_frag_atoms(struct efp *efp, size_t frag_idx, size_t n_atoms, struct efp_atom *atoms); +/** Copies information about multipole point pt_idx at fragment frag_idx into + * efp_mult_pt structure mult_pt + * + * @param[in] efp + * @param[in] frag_idx Fragment index + * @param[in] pt_idx Multipole point index + * @param[out] mult_pt efp_mult_pt to be returned + * @return ::EFP_RESULT_SUCCESS on success or error code otherwise. + */ +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. * @@ -1254,6 +1419,68 @@ void n_symm_frag(struct efp *efp, size_t *symm_frag); //static enum efp_result //check_frag_atoms(const struct frag *frag, const struct frag *lib); +/** + * Prints information of fragment atoms + * @param efp + * @param frag_index fragment index + * @param atom_index atom index in the fragment + */ +void print_atoms(struct efp *efp, size_t frag_index, size_t atom_index); + +/** + * Prints information on multipole point + * @param efp + * @param frag_index fragment index + * @param pt_index multipole point index + */ +void print_mult_pt(struct efp *efp, size_t frag_index, size_t pt_index); + +/** + * Prints information on polarizable point + * @param efp + * @param frag_index fragment index + * @param pol_index index of polarizable point + */ +void print_pol_pt(struct efp *efp, size_t frag_index, size_t pol_index); + +/** + * Prints information about ligand if any + * @param efp + * @param frag_index index of fragment + */ +void print_ligand(struct efp *efp, size_t frag_index); + +/** + * Prints information on fragment + * @param efp + * @param frag_index fragment index + */ +void print_frag_info(struct efp *efp, size_t frag_index); + +/** + * Prints information of efp_mult_pt object + * @param pt + */ +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" */ diff --git a/src/elec.c b/src/elec.c index d8a83b71..44c5d923 100644 --- a/src/elec.c +++ b/src/elec.c @@ -28,6 +28,9 @@ #include "elec.h" #include "private.h" +#include +#include + static double get_screen_damping(double r_ij, double pi, double pj) { @@ -61,106 +64,6 @@ get_screen_damping_grad(double r_ij, double pi, double pj) } } -static double -atom_mult_energy(struct efp *efp, struct frag *fr_i, struct frag *fr_j, - size_t atom_i_idx, size_t pt_j_idx, const struct swf *swf) -{ - struct efp_atom *at_i = fr_i->atoms + atom_i_idx; - struct multipole_pt *pt_j = fr_j->multipole_pts + pt_j_idx; - - vec_t dr = { - pt_j->x - at_i->x - swf->cell.x, - pt_j->y - at_i->y - swf->cell.y, - pt_j->z - at_i->z - swf->cell.z - }; - - double energy = 0.0, ccdamp = 1.0; - - if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { - double r = vec_len(&dr); - double sp = fr_j->screen_params[pt_j_idx]; - - ccdamp = get_screen_damping(r, sp, HUGE_VAL); - } - - /* charge - monopole */ - energy += ccdamp * efp_charge_charge_energy(at_i->znuc, - pt_j->monopole, &dr); - - /* charge - dipole */ - energy += efp_charge_dipole_energy(at_i->znuc, &pt_j->dipole, &dr); - - /* charge - quadrupole */ - energy += efp_charge_quadrupole_energy(at_i->znuc, - pt_j->quadrupole, &dr); - - /* charge - octupole */ - energy += efp_charge_octupole_energy(at_i->znuc, pt_j->octupole, &dr); - - return energy; -} - -static void -atom_mult_grad(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, - size_t atom_i_idx, size_t pt_j_idx, const struct swf *swf) -{ - const struct frag *fr_i = efp->frags + fr_i_idx; - const struct frag *fr_j = efp->frags + fr_j_idx; - const struct efp_atom *at_i = fr_i->atoms + atom_i_idx; - const struct multipole_pt *pt_j = fr_j->multipole_pts + pt_j_idx; - - vec_t dr = { - pt_j->x - at_i->x - swf->cell.x, - pt_j->y - at_i->y - swf->cell.y, - pt_j->z - at_i->z - swf->cell.z - }; - - vec_t force_, torque_i_, torque_j_; - vec_t force = vec_zero, torque_i = vec_zero, torque_j = vec_zero; - - /* charge - charge */ - efp_charge_charge_grad(at_i->znuc, pt_j->monopole, &dr, - &force_, &torque_i_, &torque_j_); - - if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { - double r = vec_len(&dr); - double sp = fr_j->screen_params[pt_j_idx]; - double gdamp = get_screen_damping_grad(r, sp, HUGE_VAL); - - force_.x *= gdamp; - force_.y *= gdamp; - force_.z *= gdamp; - } - - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* charge - dipole */ - efp_charge_dipole_grad(at_i->znuc, &pt_j->dipole, &dr, - &force_, &torque_i_, &torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* charge - quadrupole */ - efp_charge_quadrupole_grad(at_i->znuc, pt_j->quadrupole, &dr, - &force_, &torque_i_, &torque_j_); - vec_negate(&torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* charge - octupole */ - efp_charge_octupole_grad(at_i->znuc, pt_j->octupole, &dr, - &force_, &torque_i_, &torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - vec_scale(&force, swf->swf); - vec_scale(&torque_i, swf->swf); - vec_scale(&torque_j, swf->swf); - - efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(at_i->x), - &force, &torque_i); - efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), - &force, &torque_j); - efp_add_stress(&swf->dr, &force, &efp->stress); -} - static double mult_mult_energy(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, size_t pt_i_idx, size_t pt_j_idx, const struct swf *swf) @@ -176,56 +79,79 @@ mult_mult_energy(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, pt_j->z - pt_i->z - swf->cell.z }; - double energy = 0.0, ccdamp = 1.0; + double energy = 0.0, ccdamp = 1.0, ccdamp_i = 1.0, ccdamp_j = 1.0; + double qi = pt_i->monopole + pt_i->znuc; + double qj = pt_j->monopole + pt_j->znuc; if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { double r = vec_len(&dr); - double screen_i = fr_i->screen_params[pt_i_idx]; - double screen_j = fr_j->screen_params[pt_j_idx]; + //double screen_i = fr_i->screen_params[pt_i_idx]; + //double screen_j = fr_j->screen_params[pt_j_idx]; + double screen_i = pt_i->screen2; + double screen_j = pt_j->screen2; ccdamp = get_screen_damping(r, screen_i, screen_j); + ccdamp_j = get_screen_damping(r, screen_j, HUGE_VAL); + ccdamp_i = get_screen_damping(r, screen_i, HUGE_VAL); } - /* monopole - monopole */ - energy += ccdamp * efp_charge_charge_energy(pt_i->monopole, - pt_j->monopole, &dr); + // charge - monopole + if (pt_i->if_znuc && pt_j->if_mon) + energy += ccdamp_j * efp_charge_charge_energy(pt_i->znuc, pt_j->monopole, &dr); + + // monopole - charge + if (pt_j->if_znuc && pt_i->if_mon) + energy += ccdamp_i * efp_charge_charge_energy(pt_j->znuc, pt_i->monopole, &dr); + + // charge-charge + if (pt_i->if_znuc && pt_j->if_znuc) + energy += efp_charge_charge_energy(pt_i->znuc, pt_j->znuc, &dr); + + // monopole - monopole + if (pt_i->if_mon && pt_j->if_mon) + energy += ccdamp * efp_charge_charge_energy(pt_i->monopole, + pt_j->monopole, &dr); - /* monopole - dipole */ - energy += efp_charge_dipole_energy(pt_i->monopole, &pt_j->dipole, &dr); + // monopole - dipole + if ((pt_i->if_znuc || pt_i->if_mon) && pt_j->if_dip) + energy += efp_charge_dipole_energy(qi, &pt_j->dipole, &dr); - /* dipole - monopole */ - energy -= efp_charge_dipole_energy(pt_j->monopole, &pt_i->dipole, &dr); + // dipole - monopole + if ((pt_j->if_znuc || pt_j->if_mon) && pt_i->if_dip) + energy -= efp_charge_dipole_energy(qj, &pt_i->dipole, &dr); - /* monopole - quadrupole */ - energy += efp_charge_quadrupole_energy(pt_i->monopole, - pt_j->quadrupole, &dr); + // monopole - quadrupole + if ((pt_i->if_znuc || pt_i->if_mon) && pt_j->if_quad) + energy += efp_charge_quadrupole_energy(qi, pt_j->quadrupole, &dr); - /* quadrupole - monopole */ - energy += efp_charge_quadrupole_energy(pt_j->monopole, - pt_i->quadrupole, &dr); + // quadrupole - monopole + if ((pt_j->if_znuc || pt_j->if_mon) && pt_i->if_quad) + energy += efp_charge_quadrupole_energy(qj, pt_i->quadrupole, &dr); - /* monopole - octupole */ - energy += efp_charge_octupole_energy(pt_i->monopole, - pt_j->octupole, &dr); + // monopole - octupole + if ((pt_i->if_znuc || pt_i->if_mon) && pt_j->if_oct) + energy += efp_charge_octupole_energy(qi, pt_j->octupole, &dr); - /* octupole - monopole */ - energy -= efp_charge_octupole_energy(pt_j->monopole, - pt_i->octupole, &dr); + // octupole - monopole + if ((pt_j->if_znuc || pt_j->if_mon) && pt_i->if_oct) + energy -= efp_charge_octupole_energy(qj, pt_i->octupole, &dr); /* dipole - dipole */ - energy += efp_dipole_dipole_energy(&pt_i->dipole, &pt_j->dipole, &dr); + if (pt_i->if_dip && pt_j->if_dip) + energy += efp_dipole_dipole_energy(&pt_i->dipole, &pt_j->dipole, &dr); /* dipole - quadrupole */ - energy += efp_dipole_quadrupole_energy(&pt_i->dipole, - pt_j->quadrupole, &dr); + if (pt_i->if_dip && pt_j->if_quad) + energy += efp_dipole_quadrupole_energy(&pt_i->dipole, pt_j->quadrupole, &dr); /* quadrupole - dipole */ - energy -= efp_dipole_quadrupole_energy(&pt_j->dipole, - pt_i->quadrupole, &dr); + if (pt_j->if_dip && pt_i->if_dip) + energy -= efp_dipole_quadrupole_energy(&pt_j->dipole, pt_i->quadrupole, &dr); /* quadrupole - quadrupole */ - energy += efp_quadrupole_quadrupole_energy(pt_i->quadrupole, - pt_j->quadrupole, &dr); + if (pt_i->if_quad && pt_j->if_quad) + energy += efp_quadrupole_quadrupole_energy(pt_i->quadrupole, + pt_j->quadrupole, &dr); return energy; } @@ -245,81 +171,151 @@ mult_mult_grad(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, pt_j->z - pt_i->z - swf->cell.z }; - vec_t force_, torque_i_, torque_j_; + double qi = pt_i->monopole + pt_i->znuc; + double qj = pt_j->monopole + pt_j->znuc; + + bool if_qi = pt_i->if_mon || pt_i->if_znuc; + bool if_qj = pt_j->if_mon || pt_j->if_znuc; + + vec_t force_, torque_i_, torque_j_; vec_t force = vec_zero, torque_i = vec_zero, torque_j = vec_zero; - /* monopole - monopole */ - efp_charge_charge_grad(pt_i->monopole, pt_j->monopole, &dr, - &force_, &torque_i_, &torque_j_); + // charge-charge + if (pt_i->if_znuc && pt_j->if_znuc) { + efp_charge_charge_grad(pt_i->znuc, pt_j->znuc, &dr, + &force_, &torque_i_, &torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } - if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { - double r = vec_len(&dr); - double screen_i = fr_i->screen_params[pt_i_idx]; - double screen_j = fr_j->screen_params[pt_j_idx]; - double gdamp = get_screen_damping_grad(r, screen_i, screen_j); + // monopole - monopole + if (pt_i->if_mon && pt_j->if_mon) { + efp_charge_charge_grad(pt_i->monopole, pt_j->monopole, &dr, + &force_, &torque_i_, &torque_j_); - force_.x *= gdamp; - force_.y *= gdamp; - force_.z *= gdamp; - } + if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { + double r = vec_len(&dr); + double screen_i = pt_i->screen2; + double screen_j = pt_j->screen2; + double gdamp = get_screen_damping_grad(r, screen_i, screen_j); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* monopole - dipole */ - efp_charge_dipole_grad(pt_i->monopole, &pt_j->dipole, &dr, - &force_, &torque_i_, &torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* dipole - monopole */ - efp_charge_dipole_grad(pt_j->monopole, &pt_i->dipole, &dr, - &force_, &torque_j_, &torque_i_); - vec_negate(&force_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* monopole - quadrupole */ - efp_charge_quadrupole_grad(pt_i->monopole, pt_j->quadrupole, &dr, - &force_, &torque_i_, &torque_j_); - vec_negate(&torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* quadrupole - monopole */ - efp_charge_quadrupole_grad(pt_j->monopole, pt_i->quadrupole, &dr, - &force_, &torque_j_, &torque_i_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* monopole - octupole */ - efp_charge_octupole_grad(pt_i->monopole, pt_j->octupole, &dr, - &force_, &torque_i_, &torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); - - /* octupole - monopole */ - efp_charge_octupole_grad(pt_j->monopole, pt_i->octupole, &dr, - &force_, &torque_j_, &torque_i_); - vec_negate(&force_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + force_.x *= gdamp; + force_.y *= gdamp; + force_.z *= gdamp; + } + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // charge-monopole + if (pt_i->if_znuc && pt_j->if_mon) { + efp_charge_charge_grad(pt_i->znuc, pt_j->monopole, &dr, + &force_, &torque_i_, &torque_j_); + + if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { + double r = vec_len(&dr); + double screen_j = pt_j->screen2; + double gdamp_j = get_screen_damping_grad(r, screen_j, HUGE_VAL); + + force_.x *= gdamp_j; + force_.y *= gdamp_j; + force_.z *= gdamp_j; + //printf("MM %d %d %d %d gdamp %lf \n", + // fr_i_idx, fr_j_idx, pt_i_idx, pt_j_idx, gdamp_j); + } + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // monopole-charge + if (pt_j->if_znuc && pt_i->if_mon) { + efp_charge_charge_grad(pt_i->monopole, pt_j->znuc, &dr, + &force_, &torque_i_, &torque_j_); + + if (efp->opts.elec_damp == EFP_ELEC_DAMP_SCREEN) { + double r = vec_len(&dr); + double screen_i = pt_i->screen2; + double gdamp_i = get_screen_damping_grad(r, screen_i, HUGE_VAL); + + force_.x *= gdamp_i; + force_.y *= gdamp_i; + force_.z *= gdamp_i; + } + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // monopole - dipole + if (if_qi && pt_j->if_dip) { + efp_charge_dipole_grad(qi, &pt_j->dipole, &dr, + &force_, &torque_i_, &torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // dipole - monopole + if (if_qj && pt_i->if_dip) { + efp_charge_dipole_grad(qj, &pt_i->dipole, &dr, + &force_, &torque_j_, &torque_i_); + vec_negate(&force_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // monopole - quadrupole + if (if_qi && pt_j->if_quad) { + efp_charge_quadrupole_grad(qi, pt_j->quadrupole, &dr, + &force_, &torque_i_, &torque_j_); + vec_negate(&torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // quadrupole - monopole + if (if_qj && pt_i->if_quad) { + efp_charge_quadrupole_grad(qj, pt_i->quadrupole, &dr, + &force_, &torque_j_, &torque_i_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // monopole - octupole + if (if_qi && pt_j->if_oct) { + efp_charge_octupole_grad(qi, pt_j->octupole, &dr, + &force_, &torque_i_, &torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } + + // octupole - monopole + if (if_qj && pt_i->if_oct) { + efp_charge_octupole_grad(qj, pt_i->octupole, &dr, + &force_, &torque_j_, &torque_i_); + vec_negate(&force_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } /* dipole - dipole */ - efp_dipole_dipole_grad(&pt_i->dipole, &pt_j->dipole, &dr, - &force_, &torque_i_, &torque_j_); - vec_negate(&torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + if (pt_i->if_dip && pt_j->if_dip) { + efp_dipole_dipole_grad(&pt_i->dipole, &pt_j->dipole, &dr, + &force_, &torque_i_, &torque_j_); + vec_negate(&torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } /* dipole - quadrupole */ - efp_dipole_quadrupole_grad(&pt_i->dipole, pt_j->quadrupole, &dr, - &force_, &torque_i_, &torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + if (pt_i->if_dip && pt_j->if_quad) { + efp_dipole_quadrupole_grad(&pt_i->dipole, pt_j->quadrupole, &dr, + &force_, &torque_i_, &torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } /* quadrupole - dipole */ - efp_dipole_quadrupole_grad(&pt_j->dipole, pt_i->quadrupole, &dr, - &force_, &torque_j_, &torque_i_); - vec_negate(&force_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + if (pt_j->if_dip && pt_i->if_quad) { + efp_dipole_quadrupole_grad(&pt_j->dipole, pt_i->quadrupole, &dr, + &force_, &torque_j_, &torque_i_); + vec_negate(&force_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } /* quadrupole - quadrupole */ - efp_quadrupole_quadrupole_grad(pt_i->quadrupole, pt_j->quadrupole, - &dr, &force_, &torque_i_, &torque_j_); - vec_negate(&torque_j_); - add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + if (pt_i->if_quad && pt_j->if_quad) { + efp_quadrupole_quadrupole_grad(pt_i->quadrupole, pt_j->quadrupole, + &dr, &force_, &torque_i_, &torque_j_); + vec_negate(&torque_j_); + add_3(&force, &force_, &torque_i, &torque_i_, &torque_j, &torque_j_); + } vec_scale(&force, swf->swf); vec_scale(&torque_i, swf->swf); @@ -345,65 +341,6 @@ efp_frag_frag_elec(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx) return 0.0; } else { - /* nuclei - nuclei */ - for (size_t ii = 0; ii < fr_i->n_atoms; ii++) { - for (size_t jj = 0; jj < fr_j->n_atoms; jj++) { - struct efp_atom *at_i = fr_i->atoms + ii; - struct efp_atom *at_j = fr_j->atoms + jj; - - vec_t dr = { - at_j->x - at_i->x - swf.cell.x, - at_j->y - at_i->y - swf.cell.y, - at_j->z - at_i->z - swf.cell.z - }; - - energy += efp_charge_charge_energy(at_i->znuc, - at_j->znuc, &dr); - if (efp->do_gradient) { - vec_t force, add_i, add_j; - - efp_charge_charge_grad(at_i->znuc, at_j->znuc, - &dr, &force, &add_i, &add_j); - vec_scale(&force, swf.swf); - efp_add_force(efp->grad + fr_i_idx, - CVEC(fr_i->x), CVEC(at_i->x), &force, NULL); - efp_sub_force(efp->grad + fr_j_idx, - CVEC(fr_j->x), CVEC(at_j->x), &force, NULL); - efp_add_stress(&swf.dr, &force, &efp->stress); - } - } - } - - /* nuclei - mult points */ - for (size_t ii = 0; ii < fr_i->n_atoms; ii++) { - for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) { - energy += atom_mult_energy(efp, fr_i, fr_j, - ii, jj, &swf); - if (efp->do_gradient) { - atom_mult_grad(efp, fr_i_idx, fr_j_idx, - ii, jj, &swf); - } - } - } - - /* mult points - nuclei */ - for (size_t jj = 0; jj < fr_j->n_atoms; jj++) { - for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) { - struct swf swf2 = swf; - - vec_negate(&swf2.cell); - vec_negate(&swf2.dr); - vec_negate(&swf2.dswf); - - energy += atom_mult_energy(efp, fr_j, fr_i, - jj, ii, &swf2); - if (efp->do_gradient) { - atom_mult_grad(efp, fr_j_idx, fr_i_idx, - jj, ii, &swf2); - } - } - } - /* mult points - mult points */ for (size_t ii = 0; ii < fr_i->n_multipole_pts; ii++) { for (size_t jj = 0; jj < fr_j->n_multipole_pts; jj++) { @@ -480,52 +417,56 @@ efp_update_elec(struct frag *frag) CVEC(in->x), VEC(out->x)); /* rotate dipole */ - out->dipole = mat_vec(&frag->rotmat, &in->dipole); + if (out->if_dip) + out->dipole = mat_vec(&frag->rotmat, &in->dipole); /* rotate quadrupole */ - rotate_quadrupole(&frag->rotmat, - in->quadrupole, out->quadrupole); - - /* correction for Buckingham quadrupoles */ - double *quad = out->quadrupole; - - double qtr = quad[quad_idx(0, 0)] + - quad[quad_idx(1, 1)] + - quad[quad_idx(2, 2)]; - - quad[0] = 1.5 * quad[0] - 0.5 * qtr; - quad[1] = 1.5 * quad[1] - 0.5 * qtr; - quad[2] = 1.5 * quad[2] - 0.5 * qtr; - quad[3] = 1.5 * quad[3]; - quad[4] = 1.5 * quad[4]; - quad[5] = 1.5 * quad[5]; + if (out->if_quad) { + rotate_quadrupole(&frag->rotmat, in->quadrupole, out->quadrupole); + + /* correction for Buckingham quadrupoles */ + double *quad = out->quadrupole; + + double qtr = quad[quad_idx(0, 0)] + + quad[quad_idx(1, 1)] + + quad[quad_idx(2, 2)]; + + quad[0] = 1.5 * quad[0] - 0.5 * qtr; + quad[1] = 1.5 * quad[1] - 0.5 * qtr; + quad[2] = 1.5 * quad[2] - 0.5 * qtr; + quad[3] = 1.5 * quad[3]; + quad[4] = 1.5 * quad[4]; + quad[5] = 1.5 * quad[5]; + } /* rotate octupole */ - rotate_octupole(&frag->rotmat, in->octupole, out->octupole); - - /* correction for Buckingham octupoles */ - double *oct = out->octupole; - - double otrx = oct[oct_idx(0, 0, 0)] + - oct[oct_idx(0, 1, 1)] + - oct[oct_idx(0, 2, 2)]; - double otry = oct[oct_idx(0, 0, 1)] + - oct[oct_idx(1, 1, 1)] + - oct[oct_idx(1, 2, 2)]; - double otrz = oct[oct_idx(0, 0, 2)] + - oct[oct_idx(1, 1, 2)] + - oct[oct_idx(2, 2, 2)]; - - oct[0] = 2.5 * oct[0] - 1.5 * otrx; - oct[1] = 2.5 * oct[1] - 1.5 * otry; - oct[2] = 2.5 * oct[2] - 1.5 * otrz; - oct[3] = 2.5 * oct[3] - 0.5 * otry; - oct[4] = 2.5 * oct[4] - 0.5 * otrz; - oct[5] = 2.5 * oct[5] - 0.5 * otrx; - oct[6] = 2.5 * oct[6] - 0.5 * otrz; - oct[7] = 2.5 * oct[7] - 0.5 * otrx; - oct[8] = 2.5 * oct[8] - 0.5 * otry; - oct[9] = 2.5 * oct[9]; + if (out->if_oct) { + rotate_octupole(&frag->rotmat, in->octupole, out->octupole); + + /* correction for Buckingham octupoles */ + double *oct = out->octupole; + + double otrx = oct[oct_idx(0, 0, 0)] + + oct[oct_idx(0, 1, 1)] + + oct[oct_idx(0, 2, 2)]; + double otry = oct[oct_idx(0, 0, 1)] + + oct[oct_idx(1, 1, 1)] + + oct[oct_idx(1, 2, 2)]; + double otrz = oct[oct_idx(0, 0, 2)] + + oct[oct_idx(1, 1, 2)] + + oct[oct_idx(2, 2, 2)]; + + oct[0] = 2.5 * oct[0] - 1.5 * otrx; + oct[1] = 2.5 * oct[1] - 1.5 * otry; + oct[2] = 2.5 * oct[2] - 1.5 * otrz; + oct[3] = 2.5 * oct[3] - 0.5 * otry; + oct[4] = 2.5 * oct[4] - 0.5 * otrz; + oct[5] = 2.5 * oct[5] - 0.5 * otrx; + oct[6] = 2.5 * oct[6] - 0.5 * otrz; + oct[7] = 2.5 * oct[7] - 0.5 * otrx; + oct[8] = 2.5 * oct[8] - 0.5 * otry; + oct[9] = 2.5 * oct[9]; + } } } @@ -535,6 +476,7 @@ compute_ai_elec_frag(struct efp *efp, size_t frag_idx) struct frag *fr_i = efp->frags + frag_idx; double energy = 0.0; + /* for (size_t i = 0; i < fr_i->n_atoms; i++) { for (size_t j = 0; j < efp->n_ptc; j++) { struct efp_atom *at_i = fr_i->atoms + i; @@ -543,27 +485,34 @@ compute_ai_elec_frag(struct efp *efp, size_t frag_idx) energy += efp_charge_charge_energy(at_i->znuc, efp->ptc[j], &dr); } - } + } */ + for (size_t i = 0; i < fr_i->n_multipole_pts; i++) { for (size_t j = 0; j < efp->n_ptc; j++) { struct multipole_pt *pt_i = fr_i->multipole_pts + i; vec_t dr = vec_sub(CVEC(pt_i->x), efp->ptc_xyz + j); /* charge - monopole */ - energy += efp_charge_charge_energy(efp->ptc[j], - pt_i->monopole, &dr); + // WHY NOT USING SCREEN2 HERE??? + if (pt_i->if_mon || pt_i->if_znuc) { + double qi = pt_i->monopole + pt_i->znuc; + energy += efp_charge_charge_energy(efp->ptc[j], + qi, &dr); + } + // energy += efp_charge_charge_energy(efp->ptc[j], + // pt_i->monopole, &dr); /* charge - dipole */ - energy += efp_charge_dipole_energy(efp->ptc[j], - &pt_i->dipole, &dr); + if (pt_i->if_dip) + energy += efp_charge_dipole_energy(efp->ptc[j], &pt_i->dipole, &dr); /* charge - quadrupole */ - energy += efp_charge_quadrupole_energy(efp->ptc[j], - pt_i->quadrupole, &dr); + if (pt_i->if_quad) + energy += efp_charge_quadrupole_energy(efp->ptc[j], pt_i->quadrupole, &dr); /* charge - octupole */ - energy += efp_charge_octupole_energy(efp->ptc[j], - pt_i->octupole, &dr); + if (pt_i->if_oct) + energy += efp_charge_octupole_energy(efp->ptc[j], pt_i->octupole, &dr); } } return energy; @@ -576,7 +525,7 @@ compute_ai_elec_frag_grad(struct efp *efp, size_t frag_idx) vec_t force, add_i, add_j, force_, add_i_, add_j_; for (size_t i = 0; i < efp->n_ptc; i++) { - /* ab initio atom - fragment atoms */ + /* ab initio atom - fragment atoms for (size_t k = 0; k < fr_j->n_atoms; k++) { struct efp_atom *at_j = fr_j->atoms + k; vec_t dr = vec_sub(CVEC(at_j->x), efp->ptc_xyz + i); @@ -586,7 +535,7 @@ compute_ai_elec_frag_grad(struct efp *efp, size_t frag_idx) vec_atomic_add(efp->ptc_grad + i, &force); efp_sub_force(efp->grad + frag_idx, CVEC(fr_j->x), CVEC(at_j->x), &force, &add_j); - } + } */ /* ab initio atom - fragment multipoles */ for (size_t k = 0; k < fr_j->n_multipole_pts; k++) { @@ -598,29 +547,38 @@ compute_ai_elec_frag_grad(struct efp *efp, size_t frag_idx) add_j = vec_zero; /* monopole */ - efp_charge_charge_grad(efp->ptc[i], pt_j->monopole, &dr, - &force_, &add_i_, &add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_mon || pt_j->if_znuc) { + double qj = pt_j->monopole + pt_j->znuc; + efp_charge_charge_grad(efp->ptc[i], qj, &dr, + &force_, &add_i_, &add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* dipole */ - efp_charge_dipole_grad(efp->ptc[i], &pt_j->dipole, &dr, - &force_, &add_i_, &add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_dip) { + efp_charge_dipole_grad(efp->ptc[i], &pt_j->dipole, &dr, + &force_, &add_i_, &add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* quadrupole */ - efp_charge_quadrupole_grad(efp->ptc[i], - pt_j->quadrupole, &dr, &force_, &add_i_, &add_j_); - vec_negate(&add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_quad) { + efp_charge_quadrupole_grad(efp->ptc[i], + pt_j->quadrupole, &dr, &force_, &add_i_, &add_j_); + vec_negate(&add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* octupole */ - efp_charge_octupole_grad(efp->ptc[i], pt_j->octupole, - &dr, &force_, &add_i_, &add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_oct) { + efp_charge_octupole_grad(efp->ptc[i], pt_j->octupole, + &dr, &force_, &add_i_, &add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } vec_atomic_add(efp->ptc_grad + i, &force); efp_sub_force(efp->grad + frag_idx, CVEC(fr_j->x), diff --git a/src/parse.c b/src/parse.c index 8e5d2f64..3040cb37 100644 --- a/src/parse.c +++ b/src/parse.c @@ -28,10 +28,28 @@ #include #include #include +#include #include "stream.h" #include "private.h" +static void init_multipole_pt(struct multipole_pt *pt) { + memset(pt, 0, sizeof(*pt)); + pt->screen2 = 10.0; + pt->screen0 = 10.0; + pt->if_znuc = false; + pt->if_mon = false; + pt->if_dip = false; + pt->if_quad = false; + pt->if_oct = false; + pt->if_scr2 = false; + 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) { @@ -172,10 +190,14 @@ parse_coordinates(struct frag *frag, struct stream *stream) struct multipole_pt *last_pt = frag->multipole_pts + frag->n_multipole_pts - 1; - memset(last_pt, 0, sizeof(*last_pt)); + //memset(last_pt, 0, sizeof(*last_pt)); + init_multipole_pt(last_pt); + strcpy(last_pt->label,atom.label); last_pt->x = atom.x; last_pt->y = atom.y; last_pt->z = atom.z; + last_pt->znuc = atom.znuc; + last_pt->if_znuc = true; efp_stream_next_line(stream); } @@ -191,13 +213,30 @@ parse_monopoles(struct frag *frag, struct stream *stream) } efp_stream_next_line(stream); - + int counter = 0; for (size_t i = 0; i < frag->n_multipole_pts; i++) { - if (!skip_label(stream) || - !tok_double(stream, &frag->multipole_pts[i].monopole) || - !tok_double(stream, NULL)){ + + if (tok_stop(stream)) { + printf(" Found %d monopoles of %d expected in fragment %s \n", + counter, frag->n_multipole_pts, frag->name); + return EFP_RESULT_SUCCESS; + } + struct multipole_pt tmp_pt; + memset(&tmp_pt, 0, sizeof(tmp_pt)); + if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label) || + !tok_double(stream, &tmp_pt.monopole) || + !tok_double(stream, &tmp_pt.znuc)){ efp_log("parse_monopoles() failure"); return EFP_RESULT_SYNTAX_ERROR; + } + for (size_t j = 0; j < frag->n_multipole_pts; j++) { + if (!strcmp(tmp_pt.label, frag->multipole_pts[j].label)) { + // found a match + frag->multipole_pts[j].monopole = tmp_pt.monopole; + frag->multipole_pts[j].if_mon = true; + counter++; + break; + } } efp_stream_next_line(stream); } @@ -220,6 +259,38 @@ parse_dipoles(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); + int counter = 0; + for (size_t i = 0; i < frag->n_multipole_pts; i++) { + + if (tok_stop(stream)) { + printf(" Found %d dipoles of %d expected in fragment %s \n", + counter, frag->n_multipole_pts, frag->name); + return EFP_RESULT_SUCCESS; + } + struct multipole_pt tmp_pt; + memset(&tmp_pt, 0, sizeof(tmp_pt)); + if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label) || + !tok_double(stream, &tmp_pt.dipole.x) || + !tok_double(stream, &tmp_pt.dipole.y) || + !tok_double(stream, &tmp_pt.dipole.z)){ + efp_log("parse_dipoles() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + for (size_t j = 0; j < frag->n_multipole_pts; j++) { + if (!strcmp(tmp_pt.label, frag->multipole_pts[j].label)) { + // found a match + frag->multipole_pts[j].dipole.x = tmp_pt.dipole.x; + frag->multipole_pts[j].dipole.y = tmp_pt.dipole.y; + frag->multipole_pts[j].dipole.z = tmp_pt.dipole.z; + frag->multipole_pts[j].if_dip = true; + counter++; + break; + } + } + efp_stream_next_line(stream); + } + +/* for (size_t i = 0; i < frag->n_multipole_pts; i++) { if (!skip_label(stream) || !tok_double(stream, &frag->multipole_pts[i].dipole.x) || @@ -230,7 +301,7 @@ parse_dipoles(struct frag *frag, struct stream *stream) } efp_stream_next_line(stream); } - +*/ if (!tok_stop(stream)){ efp_log("parse_dipoles() failure"); return EFP_RESULT_SYNTAX_ERROR; @@ -249,7 +320,41 @@ parse_quadrupoles(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); - for (size_t i = 0; i < frag->n_multipole_pts; i++) { + int counter = 0; + for (size_t i = 0; i < frag->n_multipole_pts; i++) { + + if (tok_stop(stream)) { + printf(" Found %d quadrupoles of %d expected in fragment %s \n", + counter, frag->n_multipole_pts, frag->name); + return EFP_RESULT_SUCCESS; + } + struct multipole_pt tmp_pt; + memset(&tmp_pt, 0, sizeof(tmp_pt)); + if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label) || + !tok_double(stream, &tmp_pt.quadrupole[0]) || + !tok_double(stream, &tmp_pt.quadrupole[1]) || + !tok_double(stream, &tmp_pt.quadrupole[2]) || + !tok_double(stream, &tmp_pt.quadrupole[3]) || + !tok_double(stream, &tmp_pt.quadrupole[4]) || + !tok_double(stream, &tmp_pt.quadrupole[5]) ){ + efp_log("parse_quadrupoles() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + for (size_t j = 0; j < frag->n_multipole_pts; j++) { + if (!strcmp(tmp_pt.label, frag->multipole_pts[j].label)) { + // found a match + for (int q=0; q<6; q++) { + frag->multipole_pts[j].quadrupole[q] = tmp_pt.quadrupole[q]; + } + frag->multipole_pts[j].if_quad = true; + counter++; + break; + } + } + efp_stream_next_line(stream); + } +/* + for (size_t i = 0; i < frag->n_multipole_pts; i++) { if (!skip_label(stream)){ efp_log("quadrupoles() failure"); return EFP_RESULT_SYNTAX_ERROR; @@ -265,7 +370,7 @@ parse_quadrupoles(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); } - +*/ if (!tok_stop(stream)){ efp_log("quadrupoles() failure"); return EFP_RESULT_SYNTAX_ERROR; @@ -284,6 +389,42 @@ parse_octupoles(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); + int counter = 0; + for (size_t i = 0; i < frag->n_multipole_pts; i++) { + + if (tok_stop(stream)) { + printf(" Found %d octupoles of %d expected in fragment %s \n", + counter, frag->n_multipole_pts, frag->name); + return EFP_RESULT_SUCCESS; + } + struct multipole_pt tmp_pt; + memset(&tmp_pt, 0, sizeof(tmp_pt)); + + if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label)){ + efp_log("parse_octupoles() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + + for (size_t k = 0; k < 10; k++) + if (!tok_double(stream, &tmp_pt.octupole[k])){ + efp_log("parse_octupoles() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + + for (size_t j = 0; j < frag->n_multipole_pts; j++) { + if (!strcmp(tmp_pt.label, frag->multipole_pts[j].label)) { + // found a match + for (int q=0; q<10; q++) { + frag->multipole_pts[j].octupole[q] = tmp_pt.octupole[q]; + } + frag->multipole_pts[j].if_oct = true; + counter++; + break; + } + } + efp_stream_next_line(stream); + } +/* for (size_t i = 0; i < frag->n_multipole_pts; i++) { if (!skip_label(stream)){ efp_log("parse_octupoles() failure"); @@ -300,7 +441,7 @@ parse_octupoles(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); } - +*/ if (!tok_stop(stream)){ efp_log("parse_octupoles() failure"); return EFP_RESULT_SYNTAX_ERROR; @@ -327,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"); @@ -935,7 +1078,6 @@ parse_dipquad_polarizable_pts(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); } } - printf("\n in dipquad 10"); if (!tok_stop(stream)){ efp_log("parse_dipquad_polarizable_pts() failure"); return EFP_RESULT_SYNTAX_ERROR; @@ -944,6 +1086,7 @@ parse_dipquad_polarizable_pts(struct frag *frag, struct stream *stream) return EFP_RESULT_SUCCESS; } +/* static enum efp_result parse_screen(struct frag *frag, struct stream *stream) { @@ -991,6 +1134,76 @@ parse_screen(struct frag *frag, struct stream *stream) free(scr); return EFP_RESULT_SUCCESS; } +*/ + +static enum efp_result +parse_screen(struct frag *frag, struct stream *stream) +{ + if (!frag->multipole_pts){ + efp_log("parse_screen() failure: no multipole_pts"); + return EFP_RESULT_SYNTAX_ERROR; + } + + char type; + int screen_type = -1; + type = efp_stream_get_char(stream); + if (type == '\0' || isspace(type)) { + screen_type = 0; + } + else if (type == '2') { + screen_type = 2; + } + else { + printf(" Unknown SCREEN type found for fragment %s \n", frag->name); + } + + efp_stream_next_line(stream); + + int counter = 0; + for (size_t i = 0; i < frag->n_multipole_pts; i++) { + + if (tok_stop(stream)) { + printf(" Found %d SCREEN_ parameters, %d expected in fragment %s \n", + counter, frag->n_multipole_pts, frag->name); + return EFP_RESULT_SUCCESS; + } + struct multipole_pt tmp_pt; + memset(&tmp_pt, 0, sizeof(tmp_pt)); + double tmp_screen; + if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label) || + !tok_double(stream, NULL) || + !tok_double(stream, &tmp_screen)){ + efp_log("parse_screen() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + + for (size_t j = 0; j < frag->n_multipole_pts; j++) { + if (!strcmp(tmp_pt.label, frag->multipole_pts[j].label)) { + // found a match + if (screen_type == 2) { + frag->multipole_pts[j].screen2 = tmp_screen; + frag->multipole_pts[j].if_scr2 = true; + //printf(" SCREEN2 param %lf \n", frag->multipole_pts[j].screen2); + } + if (screen_type == 0) { + frag->multipole_pts[j].screen0 = tmp_screen; + frag->multipole_pts[j].if_scr0 = true; + //printf(" SCREEN0 param %lf \n", frag->multipole_pts[j].screen0); + } + counter++; + break; + } + } + efp_stream_next_line(stream); + } + + if (!tok_stop(stream)){ + efp_log("parse_screen() failure"); + return EFP_RESULT_SYNTAX_ERROR; + } + + return EFP_RESULT_SUCCESS; +} static enum efp_result parse_xrfit(struct frag *frag, struct stream *stream) @@ -1069,7 +1282,7 @@ get_parse_fn(struct stream *stream) { "CTVEC", skip_ctvec }, { "CTFOK", skip_ctfok }, { "DIPOLE-QUADRUPOLE DYNAMIC POLARIZABLE POINTS", parse_dipquad_polarizable_pts}, - { "SCREEN", parse_screen }, + { "SCREEN", parse_screen }, { "XRFIT", parse_xrfit }, { "POLAB", parse_polab }, }; @@ -1184,3 +1397,5 @@ efp_add_potential(struct efp *efp, const char *path) return res; } + + diff --git a/src/pol.c b/src/pol.c index 4dc85755..98a56cc5 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) { @@ -81,34 +75,40 @@ get_multipole_field(const vec_t *xyz, const struct multipole_pt *mult_pt, double r7 = r5 * r * r; /* charge */ - field.x += swf->swf * mult_pt->monopole * dr.x / r3; - field.y += swf->swf * mult_pt->monopole * dr.y / r3; - field.z += swf->swf * mult_pt->monopole * dr.z / r3; + if (mult_pt->if_mon || mult_pt->if_dip) { + field.x += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.x / r3; + field.y += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.y / r3; + field.z += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.z / r3; + } /* dipole */ - t1 = vec_dot(&mult_pt->dipole, &dr); + if (mult_pt->if_dip) { + t1 = vec_dot(&mult_pt->dipole, &dr); - field.x += swf->swf * (3.0 / r5 * t1 * dr.x - mult_pt->dipole.x / r3); - field.y += swf->swf * (3.0 / r5 * t1 * dr.y - mult_pt->dipole.y / r3); - field.z += swf->swf * (3.0 / r5 * t1 * dr.z - mult_pt->dipole.z / r3); + field.x += swf->swf * (3.0 / r5 * t1 * dr.x - mult_pt->dipole.x / r3); + field.y += swf->swf * (3.0 / r5 * t1 * dr.y - mult_pt->dipole.y / r3); + field.z += swf->swf * (3.0 / r5 * t1 * dr.z - mult_pt->dipole.z / r3); + } /* quadrupole */ - t1 = quadrupole_sum(mult_pt->quadrupole, &dr); - - t2 = mult_pt->quadrupole[quad_idx(0, 0)] * dr.x + - mult_pt->quadrupole[quad_idx(1, 0)] * dr.y + - mult_pt->quadrupole[quad_idx(2, 0)] * dr.z; - field.x += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.x); - - t2 = mult_pt->quadrupole[quad_idx(0, 1)] * dr.x + - mult_pt->quadrupole[quad_idx(1, 1)] * dr.y + - mult_pt->quadrupole[quad_idx(2, 1)] * dr.z; - field.y += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.y); - - t2 = mult_pt->quadrupole[quad_idx(0, 2)] * dr.x + - mult_pt->quadrupole[quad_idx(1, 2)] * dr.y + - mult_pt->quadrupole[quad_idx(2, 2)] * dr.z; - field.z += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.z); + if (mult_pt->if_quad) { + t1 = quadrupole_sum(mult_pt->quadrupole, &dr); + + t2 = mult_pt->quadrupole[quad_idx(0, 0)] * dr.x + + mult_pt->quadrupole[quad_idx(1, 0)] * dr.y + + mult_pt->quadrupole[quad_idx(2, 0)] * dr.z; + field.x += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.x); + + t2 = mult_pt->quadrupole[quad_idx(0, 1)] * dr.x + + mult_pt->quadrupole[quad_idx(1, 1)] * dr.y + + mult_pt->quadrupole[quad_idx(2, 1)] * dr.z; + field.y += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.y); + + t2 = mult_pt->quadrupole[quad_idx(0, 2)] * dr.x + + mult_pt->quadrupole[quad_idx(1, 2)] * dr.y + + mult_pt->quadrupole[quad_idx(2, 2)] * dr.z; + field.z += swf->swf * (-2.0 / r5 * t2 + 5.0 / r7 * t1 * dr.z); + } /* octupole-polarizability interactions are ignored */ @@ -133,16 +133,20 @@ get_multipole_elec_potential(const vec_t *xyz, const struct multipole_pt *mult_p double r7 = r5 * r * r; /* charge */ - elpot += swf->swf * mult_pt->monopole / r; + if (mult_pt->if_mon || mult_pt->if_dip) + elpot += swf->swf * (mult_pt->monopole + mult_pt->znuc) / r; /* dipole */ - elpot += swf->swf * vec_dot(&mult_pt->dipole, &dr) / r3; + if (mult_pt->if_dip) + elpot += swf->swf * vec_dot(&mult_pt->dipole, &dr) / r3; /* quadrupole */ - elpot += swf->swf * quadrupole_sum(mult_pt->quadrupole, &dr) / r5; + if (mult_pt->if_quad) + elpot += swf->swf * quadrupole_sum(mult_pt->quadrupole, &dr) / r5; /* octupole */ - elpot += swf->swf * octupole_sum(mult_pt->octupole, &dr) / r7; + if (mult_pt->if_oct) + elpot += swf->swf * octupole_sum(mult_pt->octupole, &dr) / r7; return elpot; } @@ -165,29 +169,6 @@ get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx) if (swf.swf == 0.0) continue; - /* field 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 = { - pt->x - at->x - swf.cell.x, - pt->y - at->y - swf.cell.y, - pt->z - at->z - swf.cell.z - }; - - double r = vec_len(&dr); - double r3 = r * r * r; - double p1 = 1.0; - - if (efp->opts.pol_damp == EFP_POL_DAMP_TT) { - p1 = efp_get_pol_damp_tt(r, fr_i->pol_damp, - fr_j->pol_damp); - } - elec_field.x += swf.swf * at->znuc * dr.x / r3 * p1; - elec_field.y += swf.swf * at->znuc * dr.y / r3 * p1; - elec_field.z += swf.swf * at->znuc * dr.z / r3 * p1; - } - /* field due to multipoles */ for (size_t j = 0; j < fr_i->n_multipole_pts; j++) { const struct multipole_pt *mult_pt = @@ -231,9 +212,10 @@ get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx) return elec_field; } -/* this function computes electric field on a polarizable point pt_idx of fragment frag_idx due to other fragment ligand_idx */ +/* this function computes electric field on a polarizable point pt_idx + * of fragment frag_idx due to other fragment ligand_idx */ static vec_t -get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, size_t ligand_idx) +get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, int ligand_idx) { const struct frag *fr_j = efp->frags + frag_idx; const struct polarizable_pt *pt = fr_j->polarizable_pts + pt_idx; @@ -249,29 +231,6 @@ get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, size_t l if (swf.swf == 0) return elec_field; - /* field 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 = { - pt->x - at->x - swf.cell.x, - pt->y - at->y - swf.cell.y, - pt->z - at->z - swf.cell.z - }; - - double r = vec_len(&dr); - double r3 = r * r * r; - double p1 = 1.0; - - if (efp->opts.pol_damp == EFP_POL_DAMP_TT) { - p1 = efp_get_pol_damp_tt(r, fr_i->pol_damp, - fr_j->pol_damp); - } - elec_field.x += swf.swf * at->znuc * dr.x / r3 * p1; - elec_field.y += swf.swf * at->znuc * dr.y / r3 * p1; - elec_field.z += swf.swf * at->znuc * dr.z / r3 * p1; - } - /* field due to multipoles */ for (size_t j = 0; j < fr_i->n_multipole_pts; j++) { const struct multipole_pt *mult_pt = @@ -314,63 +273,22 @@ get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, size_t l 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) { - vec_t *elec_field = (vec_t *)data; + (void)data; #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif for (size_t i = from; i < to; i++) { - const struct frag *frag = efp->frags + i; + // const struct frag *frag = efp->frags + i; + struct frag *frag = efp->frags + i; for (size_t j = 0; j < frag->n_polarizable_pts; j++) { - elec_field[frag->polarizable_offset + j] = - get_elec_field(efp, i, j); + struct polarizable_pt *pt = frag->polarizable_pts + j; + + pt->elec_field = get_elec_field(efp, i, j); } } } @@ -378,100 +296,58 @@ compute_elec_field_range(struct efp *efp, size_t from, size_t to, void *data) static void compute_ligand_field_range(struct efp *efp, size_t from, size_t to, void *data) { - vec_t *elec_field = (vec_t *)data; + (void)data; #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif for (size_t i = from; i < to; i++) { - const struct frag *frag = efp->frags + i; + struct frag *frag = efp->frags + i; for (size_t j = 0; j < frag->n_polarizable_pts; j++) { - elec_field[frag->polarizable_offset + j] = - get_ligand_field(efp, i, j, efp->opts.ligand); - } + struct polarizable_pt *pt = frag->polarizable_pts + j; + pt->ligand_field = get_ligand_field(efp, i, j, efp->opts.ligand); + } } } static void compute_fragment_field_range(struct efp *efp, size_t from, size_t to, void *data) { - vec_t *elec_field = (vec_t *)data; - struct frag *ligand = efp->frags + efp->opts.ligand; + (void)data; + + struct ligand *ligand = efp->ligand; #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif for (size_t i = from; i < to; i++) { - const struct frag *frag = efp->frags + i; + struct frag *frag = efp->frags + i; - for (size_t j = 0; j < ligand->n_polarizable_pts; j++) { - elec_field[frag->fragment_field_offset + j] = - get_ligand_field(efp, efp->opts.ligand, j, i); + for (size_t j = 0; j < ligand->n_ligand_pts; j++) { + struct ligand_pt *pt = ligand->ligand_pts + j; + pt->fragment_field[i] = get_ligand_field(efp, efp->ligand_index, j, i); } } } static enum efp_result -compute_elec_field(struct efp *efp) -{ - vec_t *elec_field; - vec_t *ligand_field; - vec_t *fragment_field; - int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand != -1) ? 1 : 0; +compute_elec_field(struct efp *efp) { + enum efp_result res; - elec_field = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); - efp_balance_work(efp, compute_elec_field_range, elec_field); - efp_allreduce((double *)elec_field, 3 * efp->n_polarizable_pts); + efp_balance_work(efp, compute_elec_field_range, NULL); - if (do_pairwise) { + if (efp->opts.enable_pairwise) { // this is field due to ligand on a fragment point - ligand_field = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); - efp_balance_work(efp, compute_ligand_field_range, ligand_field); - efp_allreduce((double *)ligand_field, 3 * efp->n_polarizable_pts); + // 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 - fragment_field = (vec_t *) calloc(efp->n_fragment_field_pts, sizeof(vec_t)); - efp_balance_work(efp, compute_fragment_field_range, fragment_field); - efp_allreduce((double *) fragment_field, 3 * efp->n_fragment_field_pts); + efp_balance_work(efp, compute_fragment_field_range, NULL); } } - struct frag *ligand; - if (do_pairwise && efp->opts.ligand != -1) ligand = efp->frags + efp->opts.ligand; -#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 = - elec_field[frag->polarizable_offset + j]; - frag->polarizable_pts[j].elec_field_wf = vec_zero; - if (do_pairwise) { - frag->polarizable_pts[j].ligand_field = - ligand_field[frag->polarizable_offset + j]; - } - } - if (do_pairwise && i!= efp->opts.ligand && efp->opts.ligand != -1) { - for (size_t lp = 0; lp < ligand->n_polarizable_pts; lp++) { - efp->fragment_field[frag->fragment_field_offset + lp] = - fragment_field[frag->fragment_field_offset + lp]; - } - } - - } - free(elec_field); - if (do_pairwise) { - free(ligand_field); - if (efp->opts.ligand != -1) - free(fragment_field); - } - - if (efp->opts.terms & EFP_TERM_AI_POL) - if ((res = add_electron_density_field(efp))) - return res; return EFP_RESULT_SUCCESS; } @@ -479,7 +355,7 @@ compute_elec_field(struct efp *efp) static enum efp_result compute_elec_field_crystal(struct efp *efp) { - int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand != -1) ? 1 : 0; + int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand > -1) ? 1 : 0; enum efp_result res; int nsymm = efp->nsymm_frag; @@ -491,30 +367,26 @@ compute_elec_field_crystal(struct efp *efp) struct frag *frag = efp->frags + unique_frag[i]; for (size_t j = 0; j < frag->n_polarizable_pts; j++) { - //elec_field[frag->polarizable_offset + j] = get_elec_field(efp, i, j); frag->polarizable_pts[j].elec_field = get_elec_field(efp, unique_frag[i], j); frag->polarizable_pts[j].elec_field_wf = vec_zero; } } - // slow and painful if (do_pairwise) { - - struct frag *ligand; - ligand = efp->frags + efp->opts.ligand; + struct ligand *ligand = efp->ligand; for (size_t i = 0; i < efp->n_frag; i++) { - if (i != efp->opts.ligand) { + if (i != efp->ligand_index) { // field on ligand due to other fragments "i" struct frag *frag = efp->frags + i; - for (size_t lp = 0; lp < ligand->n_polarizable_pts; lp++) { - efp->fragment_field[frag->fragment_field_offset + lp] = - get_ligand_field(efp, efp->opts.ligand, lp, i); + for (size_t lp = 0; lp < ligand->n_ligand_pts; lp++) { + struct ligand_pt *pt = ligand->ligand_pts + lp; + pt->fragment_field[i] = get_ligand_field(efp, efp->ligand_index, lp, i); } // field on fragments "i" due to ligand for (size_t p = 0; p < frag->n_polarizable_pts; p++) { frag->polarizable_pts[p].ligand_field = - get_ligand_field(efp, i, p, efp->opts.ligand); + get_ligand_field(efp, i, p, efp->ligand_index); } } } @@ -548,7 +420,6 @@ get_induced_dipole_field(struct efp *efp, size_t frag_idx, for (size_t jj = 0; jj < fr_j->n_polarizable_pts; jj++) { struct polarizable_pt *pt_j = fr_j->polarizable_pts+jj; - size_t idx = fr_j->polarizable_offset+jj; vec_t dr = { pt->x - pt_j->x + swf.cell.x, @@ -560,8 +431,8 @@ get_induced_dipole_field(struct efp *efp, size_t frag_idx, double r3 = r * r * r; double r5 = r3 * r * r; - double t1 = vec_dot(&efp->indip[idx], &dr); - double t2 = vec_dot(&efp->indipconj[idx], &dr); + double t1 = vec_dot(&pt_j->indip, &dr); + double t2 = vec_dot(&pt_j->indipconj, &dr); double p1 = 1.0; @@ -569,32 +440,97 @@ get_induced_dipole_field(struct efp *efp, size_t frag_idx, p1 = efp_get_pol_damp_tt(r, fr_i->pol_damp, fr_j->pol_damp); } - field->x -= swf.swf * p1 * (efp->indip[idx].x / r3 - + field->x -= swf.swf * p1 * (pt_j->indip.x / r3 - 3.0 * t1 * dr.x / r5); - field->y -= swf.swf * p1 * (efp->indip[idx].y / r3 - + field->y -= swf.swf * p1 * (pt_j->indip.y / r3 - 3.0 * t1 * dr.y / r5); - field->z -= swf.swf * p1 * (efp->indip[idx].z / r3 - + field->z -= swf.swf * p1 * (pt_j->indip.z / r3 - 3.0 * t1 * dr.z / r5); field_conj->x -= swf.swf * p1 * - (efp->indipconj[idx].x / r3 - 3.0 * t2 * dr.x / r5); + (pt_j->indipconj.x / r3 - 3.0 * t2 * dr.x / r5); field_conj->y -= swf.swf * p1 * - (efp->indipconj[idx].y / r3 - 3.0 * t2 * dr.y / r5); + (pt_j->indipconj.y / r3 - 3.0 * t2 * dr.y / r5); field_conj->z -= swf.swf * p1 * - (efp->indipconj[idx].z / r3 - 3.0 * t2 * dr.z / r5); - } + (pt_j->indipconj.z / r3 - 3.0 * t2 * dr.z / r5); + } } } +static bool +if_clean_indip(struct efp *efp, int counter) { + return (counter == 0); +} + +static bool +if_clean_field(struct efp *efp, int counter) { + return (counter == 0); +} + +static void +zero_ind_dipoles(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 = vec_zero; + pt->indipconj = vec_zero; + pt->indip_old = vec_zero; + pt->indipconj_old = vec_zero; + } + } +} + +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) +{ + +#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->elec_field = vec_zero; + pt->ligand_field = vec_zero; + pt->elec_field_wf = vec_zero; + } + } +} static void compute_id_range(struct efp *efp, size_t from, size_t to, void *data) { double conv = 0.0; - vec_t *id_new, *id_conj_new; - - id_new = ((struct id_work_data *)data)->id_new; - id_conj_new = ((struct id_work_data *)data)->id_conj_new; #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) reduction(+:conv) @@ -604,7 +540,6 @@ compute_id_range(struct efp *efp, size_t from, size_t to, void *data) for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; - size_t idx = frag->polarizable_offset + j; vec_t field, field_conj; /* electric field from other induced dipoles */ @@ -620,18 +555,18 @@ compute_id_range(struct efp *efp, size_t from, size_t to, void *data) field_conj.y += pt->elec_field.y + pt->elec_field_wf.y; field_conj.z += pt->elec_field.z + pt->elec_field_wf.z; - id_new[idx] = mat_vec(&pt->tensor, &field); - id_conj_new[idx] = mat_trans_vec(&pt->tensor, - &field_conj); + memcpy(&pt->indip_old, &pt->indip, sizeof(vec_t)); + memcpy(&pt->indipconj_old, &pt->indipconj, sizeof(vec_t)); + + pt->indip = mat_vec(&pt->tensor, &field); + pt->indipconj = mat_trans_vec(&pt->tensor, &field_conj); - conv += vec_dist(&id_new[idx], &efp->indip[idx]); - conv += vec_dist(&id_conj_new[idx], - &efp->indipconj[idx]); + conv += vec_dist(&pt->indip, &pt->indip_old); + conv += vec_dist(&pt->indipconj, &pt->indipconj_old); } } - ((struct id_work_data *)data)->conv += conv; - + *(double *)data += conv; } static void @@ -639,10 +574,6 @@ compute_id_crystal(struct efp *efp, void *data) { // This is not parallelized!!! double conv = 0.0; - vec_t *id_new, *id_conj_new; - - id_new = ((struct id_work_data *)data)->id_new; - id_conj_new = ((struct id_work_data *)data)->id_conj_new; int nsymm = efp->nsymm_frag; size_t *unique_frag = (size_t *)calloc(nsymm, sizeof(size_t)); @@ -657,7 +588,6 @@ compute_id_crystal(struct efp *efp, void *data) for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; - size_t idx = frag->polarizable_offset + j; vec_t field, field_conj; /* electric field from other induced dipoles */ @@ -672,38 +602,28 @@ compute_id_crystal(struct efp *efp, void *data) field_conj.y += pt->elec_field.y + pt->elec_field_wf.y; field_conj.z += pt->elec_field.z + pt->elec_field_wf.z; - id_new[idx] = mat_vec(&pt->tensor, &field); - id_conj_new[idx] = mat_trans_vec(&pt->tensor, - &field_conj); + memcpy(&pt->indip_old, &pt->indip, sizeof(vec_t)); + memcpy(&pt->indipconj_old, &pt->indipconj, sizeof(vec_t)); - conv += nsymm_frag[i]*vec_dist(&id_new[idx], &efp->indip[idx]); - conv += nsymm_frag[i]*vec_dist(&id_conj_new[idx], &efp->indipconj[idx]); - } - } - ((struct id_work_data *)data)->conv += conv; + pt->indip = mat_vec(&pt->tensor, &field); + pt->indipconj = mat_trans_vec(&pt->tensor, &field_conj); - // step 2: copy computed induced dipoles to fragment info - for (int i = 0; i < nsymm; i++) { - struct frag *frag = efp->frags + unique_frag[i]; - - for (size_t j = 0; j < frag->n_polarizable_pts; j++) { - size_t idx = frag->polarizable_offset + j; - - efp->indip[idx] = id_new[idx]; - efp->indipconj[idx] = id_conj_new[idx]; + conv += nsymm_frag[i]*vec_dist(&pt->indip, &pt->indip_old); + conv += nsymm_frag[i]*vec_dist(&pt->indipconj, &pt->indipconj_old); } } + *(double *)data += conv; - // step 3: broadcast induced dipoles to symmetry-identical fragments + // step 2: broadcast induced dipoles to symmetry-identical fragments for (int j=0; jfrags + unique_frag[j]; - size_t poloff_s = symmfrag->polarizable_offset; for (size_t i=0; in_frag; i++){ if (i == unique_frag[j]) // do nothing for self continue; + // found same-symmetry fragment, copy induced dipoles - //printf("\n symm i %d, symm j %d", efp->symmlist[i], efp->symmlist[unique_frag[j]]); + // printf("\n symm i %d, symm j %d", efp->symmlist[i], efp->symmlist[unique_frag[j]]); if (efp->symmlist[i] == efp->symmlist[unique_frag[j]]){ struct frag *frag = efp->frags + i; // compute rotation matrix between two fragments @@ -712,135 +632,110 @@ compute_id_crystal(struct efp *efp, void *data) //printf("\n rotmat: %lf, %lf, %lf ", rotmat.xy, rotmat.xz, rotmat.yz); for (size_t p=0; pn_polarizable_pts; p++) { - size_t poloff_i = frag->polarizable_offset; - + struct polarizable_pt *pt_symm = symmfrag->polarizable_pts + p; + struct polarizable_pt *pt = frag->polarizable_pts + p; // rotate induced dipole and copy - efp->indip[poloff_i+p] = mat_vec(&rotmat, &efp->indip[poloff_s+p]); - efp->indipconj[poloff_i+p] = mat_vec(&rotmat, &efp->indipconj[poloff_s+p]); + pt->indip = mat_vec(&rotmat, &pt_symm->indip); + pt->indipconj = mat_vec(&rotmat, &pt_symm->indipconj); } } } } - - /* - printf("\n induced dipoles"); - for (int i=0; in_polarizable_pts; i++) { - printf("\npoint %d: %lf ", i, vec_len(&efp->indip[i])); - } - */ } static double pol_scf_iter(struct efp *efp) { - struct id_work_data data; size_t npts = efp->n_polarizable_pts; - - data.conv = 0.0; - data.id_new = (vec_t *)calloc(npts, sizeof(vec_t)); - data.id_conj_new = (vec_t *)calloc(npts, sizeof(vec_t)); + double conv = 0.0; if (efp->opts.symmetry == 0) { // original case - efp_balance_work(efp, compute_id_range, &data); - - efp_allreduce((double *) data.id_new, 3 * npts); - efp_allreduce((double *) data.id_conj_new, 3 * npts); - efp_allreduce(&data.conv, 1); - - memcpy(efp->indip, data.id_new, npts * sizeof(vec_t)); - memcpy(efp->indipconj, data.id_conj_new, npts * sizeof(vec_t)); + efp_balance_work(efp, compute_id_range, &conv); + efp_allreduce(&conv, 1); } else { // crystal symmetry - compute_id_crystal(efp, &data); + compute_id_crystal(efp, &conv); } - free(data.id_new); - free(data.id_conj_new); - - for (size_t i=0; in_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; - size_t idx = frag->polarizable_offset + j; - if (vec_len(&efp->indip[idx]) > INDIP_PRINT_TRESH) { - printf("\n WARNING: induced dipole %d on fragment %d: %lf ", j, i, vec_len(&efp->indip[idx])); + // printing out information on convergence + if (efp->opts.print > 0) + printf(" IND DIPOLES NORM: %lf \n", conv / npts / 2); + if (efp->opts.print > 1) { + 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; + if (vec_len(&pt->indip) > INDIP_PRINT_TRESH) { + printf("\n WARNING: induced dipole %d on fragment %d: %lf ", j, i, vec_len(&pt->indip)); + } } } - } - /* - for (int i=0; in_polarizable_pts; i++) { - if vec_len(&efp->indip[i]) > INDIP_PRINT_TRESH { - printf("\npoint %d: %lf ", i, vec_len(&efp->indip[i])); - }; - } - */ + } - return data.conv / npts / 2; + return conv / npts / 2; } 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; - struct frag *ligand; - if (do_pairwise) { - ligand = efp->frags + efp->opts.ligand; - } + // 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 (if_pairwise) + ligand = efp->ligand; #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; - 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; - size_t idx = frag->polarizable_offset + j; - - energy += 0.5 * vec_dot(&efp->indipconj[idx], - &pt->elec_field_wf) - - 0.5 * vec_dot(&efp->indip[idx], - &pt->elec_field); - - energy1 += 0.5 * vec_dot(&efp->indipconj[idx], - &pt->elec_field_wf); - energy2 += -0.5 * vec_dot(&efp->indip[idx], - &pt->elec_field); - - if (efp->opts.enable_pairwise && i != efp->opts.ligand) { - efp->pair_energies[i].polarization += - - 0.5 * vec_dot(&efp->indip[idx], &pt->ligand_field); - ene2 += - 0.5 * vec_dot(&efp->indip[idx], &pt->ligand_field); - } + + energy += 0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf) - + 0.5 * vec_dot(&pt->indip, &pt->elec_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(&efp->indipconj[idx], &pt->elec_field_wf); - ene1 += 0.5 * vec_dot(&efp->indipconj[idx], &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); } } - if (efp->opts.enable_pairwise && i != efp->opts.ligand && efp->opts.ligand != -1) { - for (size_t lp = 0; lp < ligand->n_polarizable_pts; lp++) { - struct polarizable_pt *lpt = ligand->polarizable_pts + lp; - size_t idx = ligand->polarizable_offset + 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(&efp->indip[idx], &efp->fragment_field[frag->fragment_field_offset+lp]); - } - } - //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 @@ -848,42 +743,47 @@ 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; - size_t idx = frag->polarizable_offset + j; - energy += 0.5 * vec_dot(&efp->indipconj[idx], &pt->elec_field_wf) - - 0.5 * vec_dot(&efp->indip[idx], &pt->elec_field); + 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(&efp->indipconj[idx],&efp->indipconj_old[idx]); - vec_t vec2 = vec_sub(&efp->indip[idx], &efp->indip_old[idx]); + 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(&efp->indip[idx], &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(&efp->indipconj[idx], &pt->elec_field_wf); + 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; @@ -893,7 +793,7 @@ static void compute_energy_crystal(struct efp *efp, double *polenergy) { double energy = 0.0; - int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand != -1) ? 1 : 0; + int do_pairwise = (efp->opts.enable_pairwise && efp->opts.ligand > -1) ? 1 : 0; int nsymm = efp->nsymm_frag; size_t *unique_frag = (size_t *)calloc(nsymm, sizeof(size_t)); @@ -902,6 +802,10 @@ compute_energy_crystal(struct efp *efp, double *polenergy) size_t *nsymm_frag = (size_t *)calloc(nsymm, sizeof(size_t)); n_symm_frag(efp, nsymm_frag); + struct ligand *ligand; + if (do_pairwise) + ligand = efp->ligand; + for (int k=0; kfrags + i; @@ -910,25 +814,24 @@ compute_energy_crystal(struct efp *efp, double *polenergy) for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; - size_t idx = frag->polarizable_offset + j; - energy += factor * (0.5 * vec_dot(&efp->indipconj[idx], &pt->elec_field_wf) - - 0.5 * vec_dot(&efp->indip[idx], &pt->elec_field)); + energy += factor * (0.5 * vec_dot(&pt->indipconj, &pt->elec_field_wf) - + 0.5 * vec_dot(&pt->indip, &pt->elec_field)); - if (do_pairwise && i == efp->opts.ligand) { + if (do_pairwise && i == efp->ligand_index) { for (size_t fr = 0; fr < efp->n_frag; fr++) { if (fr == i) continue; struct frag *other_frag = efp->frags + fr; efp->pair_energies[fr].polarization += - -0.5 * - vec_dot(&efp->indip[idx], &efp->fragment_field[other_frag->fragment_field_offset + j]); + -0.5 * vec_dot(&pt->indip, &ligand->ligand_pts[j].fragment_field[fr]); + //vec_dot(&pt->indip, &efp->fragment_field[other_frag->fragment_field_offset + j]); } } } - if (do_pairwise && i == efp->opts.ligand) { + if (do_pairwise && i == efp->ligand_index) { for (size_t fr = 0; fr < efp->n_frag; fr++) { if (fr == i) continue; @@ -936,10 +839,10 @@ compute_energy_crystal(struct efp *efp, double *polenergy) struct frag *other_frag = efp->frags + fr; for (size_t p = 0; p < other_frag->n_polarizable_pts; p++) { struct polarizable_pt *pt = other_frag->polarizable_pts + p; - size_t idx = other_frag->polarizable_offset + p; + // size_t idx = other_frag->polarizable_offset + p; efp->pair_energies[fr].polarization += - -0.5 * vec_dot(&efp->indip[idx], &pt->ligand_field); + -0.5 * vec_dot(&pt->indip, &pt->ligand_field); } } } @@ -950,9 +853,6 @@ compute_energy_crystal(struct efp *efp, double *polenergy) static enum efp_result efp_compute_id_iterative(struct efp *efp) { - memset(efp->indip, 0, efp->n_polarizable_pts * sizeof(vec_t)); - memset(efp->indipconj, 0, efp->n_polarizable_pts * sizeof(vec_t)); - for (size_t iter = 1; iter <= POL_SCF_MAX_ITER; iter++) { if (pol_scf_iter(efp) < POL_SCF_TOL) break; @@ -969,8 +869,14 @@ efp_compute_pol_energy(struct efp *efp, double *energy) assert(energy); - if ((res = compute_elec_field(efp))) - return res; + // counter to know when to zero out induced dipoles and static field + // need to be explored further + static int counter = 0; + + // think how to skip recomputing static field in qm scf iterations + // check on efp->do_gradient breaks gtests... + if ((res = compute_elec_field(efp))) + return res; switch (efp->opts.pol_driver) { case EFP_POL_DRIVER_ITERATIVE: @@ -984,12 +890,22 @@ efp_compute_pol_energy(struct efp *efp, double *energy) if (res) return res; + if (efp->opts.print > 1) { + for (int i = 0; i < efp->n_frag; i++) { + printf("Fragment %d\n", i); + for (int j = 0; j < efp->frags[i].n_polarizable_pts; j++) { + print_pol_pt(efp, i, j); + } + } + } + *energy = 0.0; efp_balance_work(efp, compute_energy_range, energy); efp_allreduce(energy, 1); - 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)); + efp_balance_work(efp, copy_indip_gs, NULL); + + counter++; return EFP_RESULT_SUCCESS; } @@ -1001,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: @@ -1036,12 +949,26 @@ efp_compute_pol_energy_crystal(struct efp *efp, double *energy) assert(energy); + // counter to know when to zero out induced dipoles and static field + static int counter_crystal = 0; + + // think how to skip recomputing static field in qm scf iterations + // check on efp->do_gradient breaks gtests... if (res = compute_elec_field_crystal(efp)) return res; if (res = efp_compute_id_iterative(efp)) return res; + if (efp->opts.print > 1) { + for (int i = 0; i < efp->n_frag; i++) { + printf("Fragment %d\n", i); + for (int j = 0; j < efp->frags[i].n_polarizable_pts; j++) { + print_pol_pt(efp, i, j); + } + } + } + *energy = 0.0; compute_energy_crystal(efp, energy); @@ -1053,17 +980,16 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) { const struct frag *fr_i = efp->frags + frag_idx; const struct polarizable_pt *pt_i = fr_i->polarizable_pts + pt_idx; - size_t idx_i = fr_i->polarizable_offset + pt_idx; vec_t force, add_i, add_j, force_, add_i_, add_j_; double e; - vec_t dipole_i = { - 0.5 * (efp->indip[idx_i].x + efp->indipconj[idx_i].x), - 0.5 * (efp->indip[idx_i].y + efp->indipconj[idx_i].y), - 0.5 * (efp->indip[idx_i].z + efp->indipconj[idx_i].z) - }; + vec_t dipole_i = { + 0.5 * (pt_i->indip.x + pt_i->indipconj.x), + 0.5 * (pt_i->indip.y + pt_i->indipconj.y), + 0.5 * (pt_i->indip.z + pt_i->indipconj.z), + }; - for (size_t j = 0; j < efp->n_frag; j++) { + for (size_t j = 0; j < efp->n_frag; j++) { if (j == frag_idx || efp_skip_frag_pair(efp, frag_idx, j)) continue; @@ -1075,54 +1001,6 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) /* energy without switching applied */ double energy = 0.0; - /* induced dipole - nuclei */ - for (size_t k = 0; k < fr_j->n_atoms; k++) { - struct efp_atom *at_j = fr_j->atoms + k; - - vec_t dr = { - at_j->x - pt_i->x - swf.cell.x, - at_j->y - pt_i->y - swf.cell.y, - at_j->z - pt_i->z - swf.cell.z - }; - - double p1 = 1.0, p2 = 0.0; - - if (efp->opts.pol_damp == EFP_POL_DAMP_TT) { - double r = vec_len(&dr); - - p1 = efp_get_pol_damp_tt(r, fr_i->pol_damp, - fr_j->pol_damp); - p2 = efp_get_pol_damp_tt_grad(r, fr_i->pol_damp, - fr_j->pol_damp); - } - - e = -efp_charge_dipole_energy(at_j->znuc, - &dipole_i, &dr); - efp_charge_dipole_grad(at_j->znuc, &dipole_i, &dr, - &force, &add_j, &add_i); - vec_negate(&force); - - vec_scale(&force, p1); - vec_scale(&add_i, p1); - vec_scale(&add_j, p1); - - force.x += p2 * e * dr.x; - force.y += p2 * e * dr.y; - force.z += p2 * e * dr.z; - - vec_scale(&force, swf.swf); - vec_scale(&add_i, swf.swf); - vec_scale(&add_j, swf.swf); - - efp_add_force(efp->grad + frag_idx, CVEC(fr_i->x), - CVEC(pt_i->x), &force, &add_i); - efp_sub_force(efp->grad + j, CVEC(fr_j->x), - CVEC(at_j->x), &force, &add_j); - efp_add_stress(&swf.dr, &force, &efp->stress); - - energy += p1 * e; - } - /* induced dipole - multipoles */ for (size_t k = 0; k < fr_j->n_multipole_pts; k++) { struct multipole_pt *pt_j = fr_j->multipole_pts + k; @@ -1148,31 +1026,37 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) add_i = vec_zero; add_j = vec_zero; - /* induced dipole - charge */ - e = -efp_charge_dipole_energy(pt_j->monopole, - &dipole_i, &dr); - efp_charge_dipole_grad(pt_j->monopole, &dipole_i, &dr, - &force_, &add_j_, &add_i_); - vec_negate(&force_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + /* induced dipole - charge+monopole */ + if (pt_j->if_mon || pt_j->if_znuc) { + double qj = pt_j->monopole + pt_j->znuc; + e = -efp_charge_dipole_energy(qj, &dipole_i, &dr); + efp_charge_dipole_grad(qj, &dipole_i, &dr, + &force_, &add_j_, &add_i_); + vec_negate(&force_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* induced dipole - dipole */ - e += efp_dipole_dipole_energy(&dipole_i, - &pt_j->dipole, &dr); - efp_dipole_dipole_grad(&dipole_i, &pt_j->dipole, &dr, - &force_, &add_i_, &add_j_); - vec_negate(&add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_dip) { + e += efp_dipole_dipole_energy(&dipole_i, + &pt_j->dipole, &dr); + efp_dipole_dipole_grad(&dipole_i, &pt_j->dipole, &dr, + &force_, &add_i_, &add_j_); + vec_negate(&add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* induced dipole - quadrupole */ - e += efp_dipole_quadrupole_energy(&dipole_i, - pt_j->quadrupole, &dr); - efp_dipole_quadrupole_grad(&dipole_i, pt_j->quadrupole, - &dr, &force_, &add_i_, &add_j_); - add_3(&force, &force_, &add_i, &add_i_, - &add_j, &add_j_); + if (pt_j->if_dip) { + e += efp_dipole_quadrupole_energy(&dipole_i, + pt_j->quadrupole, &dr); + efp_dipole_quadrupole_grad(&dipole_i, pt_j->quadrupole, + &dr, &force_, &add_i_, &add_j_); + add_3(&force, &force_, &add_i, &add_i_, + &add_j, &add_j_); + } /* induced dipole - octupole interactions are ignored */ @@ -1200,7 +1084,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) /* induced dipole - induced dipoles */ for (size_t jj = 0; jj < fr_j->n_polarizable_pts; jj++) { struct polarizable_pt *pt_j = fr_j->polarizable_pts+jj; - size_t idx_j = fr_j->polarizable_offset+jj; + // size_t idx_j = fr_j->polarizable_offset+jj; vec_t dr = { pt_j->x - pt_i->x - swf.cell.x, @@ -1208,11 +1092,11 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) pt_j->z - pt_i->z - swf.cell.z }; - vec_t half_dipole_i = { - 0.5 * efp->indip[idx_i].x, - 0.5 * efp->indip[idx_i].y, - 0.5 * efp->indip[idx_i].z - }; + vec_t half_dipole_i = { + 0.5 * pt_i->indip.x, + 0.5 * pt_i->indip.y, + 0.5 * pt_i->indip.z, + }; double p1 = 1.0, p2 = 0.0; @@ -1225,11 +1109,8 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) fr_j->pol_damp); } - e = efp_dipole_dipole_energy(&half_dipole_i, - &efp->indipconj[idx_j], &dr); - efp_dipole_dipole_grad(&half_dipole_i, - &efp->indipconj[idx_j], &dr, &force, - &add_i, &add_j); + e = efp_dipole_dipole_energy(&half_dipole_i, &pt_j->indipconj, &dr); + efp_dipole_dipole_grad(&half_dipole_i, &pt_j->indipconj, &dr, &force, &add_i, &add_j); vec_negate(&add_j); vec_scale(&force, p1); @@ -1351,24 +1232,6 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz, if (swf.swf == 0.0) continue; - /* field 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); - double r3 = r * r * r; - - elec_field.x += swf.swf * at->znuc * dr.x / r3; - elec_field.y += swf.swf * at->znuc * dr.y / r3; - elec_field.z += swf.swf * at->znuc * dr.z / r3; - } - /* field 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; @@ -1383,7 +1246,6 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz, /* field 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, @@ -1394,15 +1256,17 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz, double r = vec_len(&dr); double r3 = r * r * r; double r5 = r3 * r * r; - double t1 = vec_dot(&efp->indip[idx], &dr); + // double t1 = vec_dot(&efp->indip[idx], &dr); + double t1 = vec_dot(&pt_i->indip, &dr); - elec_field.x -= swf.swf * (efp->indip[idx].x / r3 - - 3.0 * t1 * dr.x / r5); - elec_field.y -= swf.swf * (efp->indip[idx].y / r3 - - 3.0 * t1 * dr.y / r5); - elec_field.z -= swf.swf * (efp->indip[idx].z / r3 - - 3.0 * t1 * dr.z / r5); - } + elec_field.x -= swf.swf * (pt_i->indip.x / r3 - + 3.0 * t1 * dr.x / r5); + elec_field.y -= swf.swf * (pt_i->indip.y / r3 - + 3.0 * t1 * dr.y / r5); + elec_field.z -= swf.swf * (pt_i->indip.z / r3 - + 3.0 * t1 * dr.z / r5); + + } } if (efp->opts.terms & EFP_TERM_AI_POL) { @@ -1444,21 +1308,6 @@ efp_get_elec_potential(struct efp *efp, size_t frag_idx, const double *xyz, 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; @@ -1468,7 +1317,7 @@ efp_get_elec_potential(struct efp *efp, size_t frag_idx, const double *xyz, /* 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; + //size_t idx = fr_i->polarizable_offset + j; vec_t dr = { xyz[0] - pt_i->x - swf.cell.x, @@ -1479,7 +1328,7 @@ efp_get_elec_potential(struct efp *efp, size_t frag_idx, const double *xyz, 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; + elpot += swf.swf * 0.5 * (vec_dot(&pt_i->indip, &dr) + vec_dot(&pt_i->indipconj, &dr)) / r3; } } diff --git a/src/poldirect.c b/src/poldirect.c index 2380a2ab..2b7f1460 100644 --- a/src/poldirect.c +++ b/src/poldirect.c @@ -31,6 +31,9 @@ double efp_get_pol_damp_tt(double, double, double); enum efp_result efp_compute_id_direct(struct efp *); +enum efp_result efp_get_induced_dipole_values(struct efp *efp, double *dip); +enum efp_result efp_get_induced_dipole_conj_values(struct efp *efp, double *dip); +enum efp_result efp_set_induced_dipole_values(struct efp *efp, double *dip, int if_conjug); static void copy_matrix(double *dst, size_t n, size_t off_i, size_t off_j, const mat_t *m) @@ -129,21 +132,23 @@ compute_lhs(const struct efp *efp, double *c, int conj) } static void -compute_rhs(const struct efp *efp, vec_t *id, int conj) +compute_rhs(struct efp *efp, int conj) { - for (size_t i = 0, idx = 0; i < efp->n_frag; i++) { - const struct frag *frag = efp->frags + i; + 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++, idx++) { - const struct polarizable_pt *pt = + for (size_t j = 0; j < frag->n_polarizable_pts; j++) { + struct polarizable_pt *pt = frag->polarizable_pts + j; vec_t field = vec_add(&pt->elec_field, &pt->elec_field_wf); if (conj) - id[idx] = mat_trans_vec(&pt->tensor, &field); + pt->indipconj = mat_trans_vec(&pt->tensor, &field); + //id[idx] = mat_trans_vec(&pt->tensor, &field); else - id[idx] = mat_vec(&pt->tensor, &field); + pt->indip = mat_vec(&pt->tensor, &field); + //id[idx] = mat_vec(&pt->tensor, &field); } } } @@ -167,30 +172,63 @@ efp_compute_id_direct(struct efp *efp) /* induced dipoles */ compute_lhs(efp, c, 0); - compute_rhs(efp, efp->indip, 0); + compute_rhs(efp, 0); transpose_matrix(c, n); + double *dip; + dip = (double *)calloc(n, sizeof(double)); + if (efp_get_induced_dipole_values(efp, dip)) { + efp_log("efp_compute_id_direct() failure"); + } + + if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv, + (double *)dip, (fortranint_t)n) != 0) { + efp_log("dgesv: error solving for induced dipoles"); + res = EFP_RESULT_FATAL; + goto error; + } + if (efp_set_induced_dipole_values(efp, dip, 0)) { + efp_log("efp_compute_id_direct() failure"); + } + +/* if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv, (double *)efp->indip, (fortranint_t)n) != 0) { efp_log("dgesv: error solving for induced dipoles"); res = EFP_RESULT_FATAL; goto error; - } + }*/ /* conjugate induced dipoles */ compute_lhs(efp, c, 1); - compute_rhs(efp, efp->indipconj, 1); + compute_rhs(efp, 1); transpose_matrix(c, n); - if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv, + if (efp_get_induced_dipole_conj_values(efp, dip)) { + efp_log("efp_compute_id_direct() failure"); + } + + if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv, + (double *)dip, (fortranint_t)n) != 0) { + efp_log("dgesv: error solving for induced dipoles"); + res = EFP_RESULT_FATAL; + goto error; + } + if (efp_set_induced_dipole_values(efp, dip, 1)) { + efp_log("efp_compute_id_direct() failure"); + } + +/* + if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv, (double *)efp->indipconj, (fortranint_t)n) != 0) { efp_log("dgesv: error solving for conjugate induced dipoles"); res = EFP_RESULT_FATAL; goto error; - } + } */ res = EFP_RESULT_SUCCESS; error: free(c); free(ipiv); + free(dip); return res; } diff --git a/src/private.h b/src/private.h index 90237fb1..c848c6ee 100644 --- a/src/private.h +++ b/src/private.h @@ -28,6 +28,7 @@ #define LIBEFP_PRIVATE_H #include +#include #include "efp.h" #include "int.h" @@ -41,19 +42,45 @@ #define ARRAY_SIZE(arr) (sizeof(arr)/sizeof(arr[0])) struct multipole_pt { - double x, y, z; + char label[32]; /**< Multipole point label. */ + double x, y, z; + double znuc; double monopole; vec_t dipole; double quadrupole[6]; double octupole[10]; + double screen2; /***/ + double screen0; /***/ + bool if_znuc; + bool if_mon; + bool if_dip; + bool if_quad; + bool if_oct; + bool if_scr2; + bool if_scr0; }; struct polarizable_pt { double x, y, z; mat_t tensor; vec_t elec_field; - vec_t elec_field_wf; - vec_t ligand_field; + 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; // 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 */ +struct ligand_pt { + /* fields of all (n_frag) fragments on this point */ + vec_t *fragment_field; + + /* number of fragments producing fields */ + size_t n_frag; }; struct dynamic_polarizable_pt { @@ -107,12 +134,6 @@ struct frag { /* number of distributed multipole points */ size_t n_multipole_pts; - /* electrostatic screening parameters */ - double *screen_params; - - /* ab initio electrostatic screening parameters */ - double *ai_screen_params; - /* polarization damping parameter */ double pol_damp; @@ -122,7 +143,7 @@ struct frag { /* number of distributed polarizability points */ size_t n_polarizable_pts; - /* dynamic polarizability points */ + /* dynamic polarizability points */ struct dynamic_polarizable_pt *dynamic_polarizable_pts; /* dipole-quadrupole dynamic polarizability points */ @@ -167,8 +188,18 @@ struct frag { /* offset of polarizable points for this fragment */ size_t polarizable_offset; - /* offset of fragment field points for this fragment */ - size_t fragment_field_offset; +}; + +/* structure derived from struct frag for describing ligand */ +struct ligand { + /* fragment */ + struct frag *ligand_frag; + + /* array of ligand points */ + struct ligand_pt *ligand_pts; + + /* number of ligand points */ + size_t n_ligand_pts; }; struct efp { @@ -190,11 +221,11 @@ struct efp { /* array with the library of fragment updated (shifted) parameters */ struct frag **lib_current; - /* callback which computes electric field from electrons */ - efp_electron_density_field_fn get_electron_density_field; + /* pointer to ligand fragment */ + struct ligand *ligand; - /* user data for get_electron_density_field */ - void *get_electron_density_field_user_data; + /* ligand index in fragment list */ + size_t ligand_index; /* user parameters for this EFP computation */ struct efp_opts opts; @@ -223,27 +254,9 @@ struct efp { /* gradient on point charges */ vec_t *ptc_grad; - /* polarization induced dipoles */ - vec_t *indip; - - /* polarization conjugate induced dipoles */ - vec_t *indipconj; - - /* polarization induced dipoles */ - vec_t *indip_old; - - /* polarization conjugate induced dipoles */ - vec_t *indipconj_old; - /* total number of polarizable points */ size_t n_polarizable_pts; - /* total number of points to store electric field due to ligand */ - size_t n_fragment_field_pts; - - /* electric field on ligand points due to fragments */ - vec_t *fragment_field; - /* number of core orbitals in ab initio subsystem */ size_t n_ai_core; diff --git a/tests/pairwise_0.in b/tests/pairwise_0.in index 31b6dc7e..ee1f41e0 100644 --- a/tests/pairwise_0.in +++ b/tests/pairwise_0.in @@ -4,6 +4,8 @@ terms elec pol disp xr disp_damp overlap fraglib_path ../fraglib enable_pairwise true +ligand 0 + fragment h2o_l 0.000 0.000 0.000 0.000 0.000 0.000 fragment h2o_l diff --git a/tests/symm_2pw_printout.in b/tests/symm_2pw_printout.in new file mode 100644 index 00000000..218b4d0c --- /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 From 2e5fc33e5e95fabde400642c56533cbab5ad199f Mon Sep 17 00:00:00 2001 From: "Slipchenko, Lyudmila V" Date: Sun, 30 Oct 2022 18:49:57 -0400 Subject: [PATCH 2/2] version update --- .DS_Store | Bin 8196 -> 10244 bytes Doxyfile | 4 ++-- src/efp.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.DS_Store b/.DS_Store index 6d4ee305730cc7aaf9843c0bee591ed3fd76eff4..654980235b4ae30e381aa2bc42b17785ce3ced82 100644 GIT binary patch delta 1575 zcmc(fUu;uV9LIm(c7N`@>)2`O#s%7~t~6{?M!VGw-G*H&7zEsq2_xu6yI#kY_O{z? zJ4&ft3^7Jfv)q6Y#Rn55W=~_9sELV6O!PrbbPwo@@xd3PJoqvrFP`=`aWP`T6DPU% zckb`^yXVjE`#I-)`23hmM2O@PK_U9W?P(XTkN^7wc-H(CYm9IG@7DQlwF>^ z+_<>VywWJt3Eha2d@0z*hxJfHOICT;+dqXX9IzouvqKr z8yW>sYLfRoD(#<8r_VvfA3%#UWWYFg8e zDswuw-f)1=%$vAmYu(b}?b_KpFf_8-yrG3j4#{}OzP%h3$F*2WS!Ekqn%g$X@n{LS`Zx>li7}xoF z%^~iW{36rE@zbhS`*oKT6j(w$_GB^^SCZy%A*o+r2gGA}bM>b4L1|F1!!J%ur&B3K z%j%|kj~o`*m%?Oa@`EuwJDO26K7VB!5tl3*A3L^|MFvj7X6mF64bd1K=km_eMY=?n z>1|q~75a#7&^PogeMjHZP5POBrQhgxx@IF4l=eS{NV&Obi+ha9b zGx&mc?nol1rADG<^1oEl&NkK+WM#v~W^6rM#ADQKAIvKDa;CA@$a z@e(d@eJ|lPT)`W76Ibi-E|!d2tzO?9#mgU=(~@azhVKF^Hg>yldEryot>&gYP2BFS w+*!-BYgVJhb=YY!Z?fvHIW5)0Y-C(7mw7$!Wv$EW4${4-b^X=5uJwob8#>us$p8QV delta 120 zcmZn(XmOBWU|?W$DortDU;r^WfEYvza8E20o2aKK$^w!H@)?Q?l5+BsfV|U51se;e zu}^H^-OSFx!okS0*+=X*yC5@A3J4^)frKl_fQ^OUnJ4qB1abhigN$XE9M3a#v#m5c FGXMf572E&- diff --git a/Doxyfile b/Doxyfile index daaa01fa..bf369ed7 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.h b/src/efp.h index 5c55c156..43244254 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 {