Skip to content

Commit 776cbf5

Browse files
authored
Merge pull request #13 from libefp2/pol_clean
Pol clean
2 parents e52ae11 + 1eef502 commit 776cbf5

20 files changed

Lines changed: 69871 additions & 44639 deletions

EFP2.C

Lines changed: 2320 additions & 0 deletions
Large diffs are not rendered by default.

EFP2.h

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
#ifndef EFP2_H_INCLUDED
2+
#define EFP2_H_INCLUDED
3+
4+
#include <vector>
5+
6+
class EFP2_impl;
7+
class InputSection;
8+
9+
class EFP2 {
10+
public:
11+
/* Initialize EFP2 instance using user's input. */
12+
void init(InputSection &);
13+
14+
/* Returns true if the current EFP2 instance was initialized. */
15+
bool initialized();
16+
17+
/* Returns if_multipole_field value. */
18+
bool get_if_multipole_field();
19+
20+
/* Returns updated if_multipole_field value. */
21+
void update_multipole_field();
22+
23+
/* Returns if_pol_field value. */
24+
bool get_if_pol_field();
25+
26+
/* Returns updated if_pol_field value. */
27+
void update_pol_field();
28+
29+
/* Reorient EFP fragments according to Q-Chem's orientation of
30+
* ab initio subsystem. */
31+
void reorient_geometry();
32+
33+
/* Print geometry of EFP subsystem. */
34+
void print_geometry();
35+
36+
/* Setup all necessary stuff for AI/EFP dispersion computation
37+
* if enabled */
38+
void setup_aiefp_dispersion();
39+
40+
/* Compute all EFP and the rest of QM/EFP interactions after
41+
* the QM SCF has converged.
42+
* \p do_grad - specifies whether to compute the gradient.
43+
* If \p do_grad is not zero than the gradient will also be
44+
* computed. */
45+
void compute(int do_grad);
46+
47+
/* Update quantum mechanical Hamiltonian with the contributions
48+
* from EFP multipoles and atoms using AOints.
49+
* \p h - a Hamiltonian matrix.
50+
* \p code - Q-Chem's ShlPrs code. */
51+
void update_mult_ints(double *h, INTEGER code);
52+
53+
/* Update quantum mechanical Hamiltonian with the contributions
54+
* from EFP induced diples using AOints.
55+
* \p h - a Hamiltonian matrix.
56+
* \p code - Q-Chem's ShlPrs code. */
57+
void update_pol_ints(double *h, INTEGER code);
58+
59+
/* Update quantum mechanical Hamiltonian with the contributions
60+
* from EFPs using AOints.
61+
* \p h - a Hamiltonian matrix.
62+
* \p code - Q-Chem's ShlPrs code. */
63+
// void update_wf(double *h, INTEGER code);
64+
65+
/* Update quantum mechanical Hamiltonian with the contributions
66+
* from EFPs using libqints.
67+
* \p h - a Hamiltonian matrix.
68+
* \p code - Q-Chem's ShlPrs code. */
69+
void update_wf_qints(double *h, INTEGER code);
70+
71+
/* Update the positions of quantum nuclei for libefp. */
72+
void update_qm_atoms();
73+
74+
/* Returns the wavefunction dependent energy.
75+
* \p w - is a density matrix.
76+
* \p n - is a total number of elements in array \p w. */
77+
double get_wf_dependent_energy(double *w, double n);
78+
79+
/* Returns the total EFP energy. */
80+
double get_total_energy();
81+
82+
/* Returns EFP correction to the energy of the excited state.
83+
* \p w - is the density matrix of the excited state WF.
84+
* \p n - is the total number of elements in array \p w.
85+
* \p Ecis - excitation energy */
86+
double get_excited_state_energy_correction(double *w, size_t n, double Ecis);
87+
88+
/* Returns gradient on quantum nuclei from EFP electrostatics.
89+
* Upon return the \p a vector will contain 3*N elements of
90+
* gradient, where N is the number of QM atoms. */
91+
void get_qm_gradient(std::vector<double> &a);
92+
93+
/* Returns current EFP gradient.
94+
* The vector \p a will contain 6*N gradient values for
95+
* fragment translation and rotation. N is the total number of
96+
* EFP fragments. */
97+
void get_gradient(std::vector<double> &a);
98+
99+
/* computes pairwise energies between QM region and EFP fragments.
100+
* \p Escf - reference AI energy (HF or excitation energy)
101+
* \p if_excited = 1/0 a switch to distinguish between the ground state (0) or excited state (1). */
102+
void get_pairwise_energy(double Estate, int if_excited);
103+
104+
/* Print EFP energy terms to stdout. */
105+
void print_energy();
106+
107+
/* Print pairwise energy decomposition to stdout.
108+
* \p if_excited = 1/0 a switch to distinguish between the ground state (0) or excited state (1). */
109+
void print_pairwise_energy(int if_excited);
110+
111+
/* Returns a EFP2 singleton instance. */
112+
static EFP2& instance() {
113+
static EFP2 instance;
114+
return instance;
115+
}
116+
117+
private:
118+
EFP2();
119+
EFP2(const EFP2 &);
120+
~EFP2();
121+
122+
EFP2_impl *impl_;
123+
};
124+
125+
#endif /* EFP2_H_INCLUDED */

efpmd/src/main.c

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -152,7 +152,7 @@ static struct cfg *make_cfg(void)
152152
cfg_add_double(cfg, "thermostat_tau", 1.0e3);
153153
cfg_add_double(cfg, "barostat_tau", 1.0e4);
154154

155-
cfg_add_int(cfg, "ligand", 0);
155+
cfg_add_int(cfg, "ligand", -100);
156156
cfg_add_bool(cfg, "enable_pairwise", false);
157157
cfg_add_bool(cfg, "print_pbc", false);
158158
cfg_add_bool(cfg, "symmetry", false);
@@ -166,6 +166,8 @@ static struct cfg *make_cfg(void)
166166
cfg_add_int(cfg, "update_params", 0);
167167
cfg_add_double(cfg, "update_params_cutoff", 0.0);
168168

169+
cfg_add_int(cfg, "print", 0);
170+
169171
return cfg;
170172
}
171173

@@ -286,12 +288,13 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
286288
.symmetry = cfg_get_bool(cfg, "symmetry"),
287289
.symm_frag = cfg_get_enum(cfg, "symm_frag"),
288290
.update_params = cfg_get_int(cfg, "update_params"),
289-
.update_params_cutoff = cfg_get_double(cfg, "update_params_cutoff")
291+
.update_params_cutoff = cfg_get_double(cfg, "update_params_cutoff"),
292+
.print = cfg_get_int(cfg, "print")
290293
};
291294

292295
if (opts.xr_cutoff == 0.0) {
293296
opts.xr_cutoff = opts.swf_cutoff;
294-
printf("xr_cutoff is set to %lf \n\n\n", opts.xr_cutoff*0.52917721092);
297+
printf("xr_cutoff is set to %lf \n\n", opts.xr_cutoff*0.52917721092);
295298
}
296299

297300
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)
306309
add_potentials(efp, cfg, sys);
307310

308311
for (size_t i = 0; i < sys->n_frags; i++)
309-
check_fail(efp_add_fragment(efp, sys->frags[i].name));
312+
check_fail(efp_add_fragment(efp, sys->frags[i].name));
310313

311314
if (sys->n_charges > 0) {
312315
double q[sys->n_charges];
@@ -331,7 +334,10 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
331334
if (cfg_get_bool(cfg, "enable_ff"))
332335
opts.terms &= ~(EFP_TERM_ELEC | EFP_TERM_POL | EFP_TERM_DISP | EFP_TERM_XR);
333336

334-
check_fail(efp_set_opts(efp, &opts));
337+
if (opts.enable_pairwise)
338+
check_fail(efp_add_ligand(efp, opts.ligand));
339+
340+
check_fail(efp_set_opts(efp, &opts));
335341
check_fail(efp_prepare(efp));
336342
check_fail(efp_set_symmlist(efp));
337343

efpmd/src/parse.c

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -253,6 +253,10 @@ struct sys *parse_input(struct cfg *cfg, const char *path)
253253
if (cfg_get_bool(cfg, "enable_pbc"))
254254
cfg_set_bool(cfg, "enable_cutoff", true);
255255

256+
// turn off pairwise calculations if ligand is not set
257+
if (cfg_get_int(cfg, "ligand") == -100 && cfg_get_bool(cfg, "enable_pairwise") == true)
258+
error("Specify ligand for pairwise calculations");
259+
256260
check_cfg(cfg);
257261
efp_stream_close(stream);
258262

fraglib/h2o_reduced.efp

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
RUNTYP=MAKEFP EFFECTIVE FRAGMENT POTENTIAL DATA FOLLOWS...
2+
FRAGNAMEEFP GENERATED AT Tue Mar 27 16:14:27 2012
3+
$H2O_REDUCED_L
4+
EFP DATA FOR FRAGNAME SCFTYP=RHF ... GENERATED WITH BASIS SET=6-311++G(3df,2p)
5+
COORDINATES (BOHR)
6+
A01O1 0.0000000000 0.1191094785 0.0000000000 15.9949100 8.0
7+
A02H2 -1.4223059670 -0.9451766865 0.0000000000 1.0078250 1.0
8+
A03H3 1.4223059670 -0.9451766865 0.0000000000 1.0078250 1.0
9+
BO21 -0.7111529835 -0.4130336040 0.0000000000 0.0000000 0.0
10+
BO31 0.7111529835 -0.4130336040 0.0000000000 0.0000000 0.0
11+
STOP
12+
MONOPOLES
13+
A01O1 -8.4173680991 8.00000
14+
A02H2 -0.3512187169 1.00000
15+
A03H3 -0.3512187169 1.00000
16+
BO21 -0.4400972335 0.00000
17+
BO31 -0.4400972335 0.00000
18+
STOP
19+
DIPOLES
20+
A01O1 0.0000000000 -0.4127062480 0.0000000000
21+
A02H2 -0.0027268078 -0.0014022974 0.0000000000
22+
A03H3 0.0027268078 -0.0014022974 0.0000000000
23+
BO21 0.2657818726 0.2358863756 0.0000000000
24+
BO31 -0.2657818726 0.2358863756 0.0000000000
25+
STOP
26+
POLARIZABLE POINTS
27+
CT1 -0.7491370634 -0.4890369778 0.0000000020
28+
2.6629294234 1.7670590537 0.7812922107 1.2832534969 >
29+
-0.0000000037 -0.0000000033 1.1482095947 -0.0000000144 >
30+
-0.0000000132
31+
CT2 0.7491370628 -0.4890369781 0.0000000040
32+
2.6629294254 1.7670590567 0.7812922113 -1.2832535009 >
33+
0.0000000076 -0.0000000067 -1.1482095971 0.0000000195 >
34+
-0.0000000177
35+
CT3 0.0000000011 0.3892153728 -0.4957986131
36+
1.4502428775 1.8574126776 2.7113709049 0.0000000046 >
37+
-0.0000000029 -0.9972382204 0.0000000048 -0.0000000034 >
38+
-1.0913249895
39+
CT4 -0.0000000005 0.3892153799 0.4957986071
40+
1.4502428978 1.8574127185 2.7113708742 -0.0000000023 >
41+
-0.0000000009 0.9972382303 -0.0000000024 -0.0000000017 >
42+
1.0913250020
43+
STOP
44+
$END

0 commit comments

Comments
 (0)