Skip to content

Commit c0d5f29

Browse files
slipchenkoslipchenko
authored
Merge pull request #9 from libefp2/elpol
Elec potential
2 parents a215c9b + 9502441 commit c0d5f29

49 files changed

Lines changed: 266 additions & 37 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎efpmd/src/common.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@ enum run_type {
5959
RUN_TYPE_OPT,
6060
RUN_TYPE_MD,
6161
RUN_TYPE_EFIELD,
62+
RUN_TYPE_ELPOT,
6263
RUN_TYPE_GTEST,
6364
RUN_TYPE_ETEST
6465
};

‎efpmd/src/efield.c‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include "common.h"
2828

2929
void sim_efield(struct state *state);
30+
void sim_elpot(struct state *state);
3031

3132
static void
3233
print_field(size_t frag_idx, const struct efp_atom *atom, const double *field)
@@ -78,3 +79,54 @@ sim_efield(struct state *state)
7879

7980
msg("ELECTRIC FIELD JOB COMPLETED SUCCESSFULLY\n");
8081
}
82+
83+
84+
static void
85+
print_elpot(const struct efp_atom *atom, const double elpot)
86+
{
87+
msg(" %s %12.8lf %12.8lf %12.8lf %12.8lf\n",
88+
atom->label,
89+
BOHR_RADIUS * atom->x,
90+
BOHR_RADIUS * atom->y,
91+
BOHR_RADIUS * atom->z,
92+
elpot);
93+
// msg(" ELEC POTENTIAL %12.8lf\n\n", elpot);
94+
}
95+
96+
void
97+
sim_elpot(struct state *state)
98+
{
99+
size_t n_frags;
100+
101+
msg("ELECTROSTATIC POTENTIAL JOB\n\n\n");
102+
103+
print_geometry(state->efp);
104+
compute_energy(state, false);
105+
print_energy(state);
106+
check_fail(efp_get_frag_count(state->efp, &n_frags));
107+
108+
msg("COORDINATES IN ANGSTROMS, ELECTROSTATIC POTENTIAL IN ATOMIC UNITS\n");
109+
msg(" ATOM X Y Z ELPOT \n\n");
110+
111+
for (size_t i = 0; i < n_frags; i++) {
112+
double elpot;
113+
struct efp_atom *atoms;
114+
size_t n_atoms;
115+
116+
check_fail(efp_get_frag_atom_count(state->efp, i, &n_atoms));
117+
atoms = xmalloc(n_atoms * sizeof(struct efp_atom));
118+
check_fail(efp_get_frag_atoms(state->efp, i, n_atoms, atoms));
119+
120+
msg("ELECTROSTATIC POTENTIAL ON FRAGMENT %zu\n", i);
121+
122+
for (size_t j = 0; j < n_atoms; j++) {
123+
check_fail(efp_get_elec_potential(state->efp, i, &atoms[j].x, &elpot));
124+
print_elpot(atoms + j, elpot);
125+
}
126+
msg("\n");
127+
128+
free(atoms);
129+
}
130+
131+
msg("ELECTROSTATIC POTENTIAL JOB COMPLETED SUCCESSFULLY\n");
132+
}

‎efpmd/src/main.c‎

Lines changed: 20 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@ void sim_hess(struct state *);
3636
void sim_opt(struct state *);
3737
void sim_md(struct state *);
3838
void sim_efield(struct state *);
39+
void sim_elpot(struct state *);
3940
void sim_gtest(struct state *);
4041
void sim_etest(struct state *);
4142

@@ -56,6 +57,7 @@ static struct cfg *make_cfg(void)
5657
"opt\n"
5758
"md\n"
5859
"efield\n"
60+
"elpot\n"
5961
"gtest\n"
6062
"etest\n",
6163
(int []) { RUN_TYPE_SP,
@@ -64,10 +66,11 @@ static struct cfg *make_cfg(void)
6466
RUN_TYPE_OPT,
6567
RUN_TYPE_MD,
6668
RUN_TYPE_EFIELD,
69+
RUN_TYPE_ELPOT,
6770
RUN_TYPE_GTEST,
6871
RUN_TYPE_ETEST});
6972

70-
cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
73+
cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_POINTS,
7174
"xyzabc\n"
7275
"points\n"
7376
"rotmat\n",
@@ -160,20 +163,22 @@ static struct cfg *make_cfg(void)
160163
static sim_fn_t get_sim_fn(enum run_type run_type)
161164
{
162165
switch (run_type) {
163-
case RUN_TYPE_SP:
164-
return sim_sp;
165-
case RUN_TYPE_GRAD:
166-
return sim_grad;
167-
case RUN_TYPE_HESS:
168-
return sim_hess;
169-
case RUN_TYPE_OPT:
170-
return sim_opt;
171-
case RUN_TYPE_MD:
172-
return sim_md;
173-
case RUN_TYPE_EFIELD:
174-
return sim_efield;
175-
case RUN_TYPE_GTEST:
176-
return sim_gtest;
166+
case RUN_TYPE_SP:
167+
return sim_sp;
168+
case RUN_TYPE_GRAD:
169+
return sim_grad;
170+
case RUN_TYPE_HESS:
171+
return sim_hess;
172+
case RUN_TYPE_OPT:
173+
return sim_opt;
174+
case RUN_TYPE_MD:
175+
return sim_md;
176+
case RUN_TYPE_EFIELD:
177+
return sim_efield;
178+
case RUN_TYPE_ELPOT:
179+
return sim_elpot;
180+
case RUN_TYPE_GTEST:
181+
return sim_gtest;
177182
case RUN_TYPE_ETEST:
178183
return sim_etest;
179184
}

‎src/efp.h‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1062,6 +1062,20 @@ enum efp_result efp_get_frag_atoms(struct efp *efp, size_t frag_idx,
10621062
enum efp_result efp_get_electric_field(struct efp *efp, size_t frag_idx,
10631063
const double *xyz, double *field);
10641064

1065+
/**
1066+
* Get electrostatic potential for a point on a fragment.
1067+
*
1068+
* @param[in] efp The efp structure.
1069+
* @param[in] frag_idx Index of a fragment. Must be a value between zero and
1070+
* the total number of fragments minus one.
1071+
* @param[in] xyz Coordinates of a point for which electric field should be
1072+
* computed.
1073+
* @param[out] elec_potential Electrostatic potential in atomic units.
1074+
* @return ::EFP_RESULT_SUCCESS on success or error code otherwise.
1075+
*/
1076+
enum efp_result efp_get_elec_potential(struct efp *efp, size_t frag_idx,
1077+
const double *xyz, double *elec_potential);
1078+
10651079
/**
10661080
* Convert rigid body torque to derivatives of energy by Euler angles.
10671081
*

‎src/elec.h‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,28 @@ quadrupole_sum(const double *quad, const vec_t *dr)
9393
return sum;
9494
}
9595

96+
static inline double
97+
octupole_sum(const double *oct, const vec_t *dr)
98+
{
99+
/* order in which octupoles are stored */
100+
enum { xxx = 0, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz };
101+
102+
double sum = 0.0;
103+
104+
sum += oct[xxx] * dr->x * dr->x * dr->x;
105+
sum += oct[yyy] * dr->y * dr->y * dr->y;
106+
sum += oct[zzz] * dr->z * dr->z * dr->z;
107+
sum += oct[xxy] * dr->x * dr->x * dr->y * 3.0;
108+
sum += oct[xxz] * dr->x * dr->x * dr->z * 3.0;
109+
sum += oct[xyy] * dr->x * dr->y * dr->y * 3.0;
110+
sum += oct[yyz] * dr->y * dr->y * dr->z * 3.0;
111+
sum += oct[xzz] * dr->x * dr->z * dr->z * 3.0;
112+
sum += oct[yzz] * dr->y * dr->z * dr->z * 3.0;
113+
sum += oct[xyz] * dr->x * dr->y * dr->z * 6.0;
114+
115+
return sum;
116+
}
117+
96118
double efp_charge_charge_energy(double, double, const vec_t *);
97119
double efp_charge_dipole_energy(double, const vec_t *, const vec_t *);
98120
double efp_charge_quadrupole_energy(double, const double *, const vec_t *);

‎src/electerms.c‎

Lines changed: 0 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -26,28 +26,6 @@
2626

2727
#include "elec.h"
2828

29-
static double
30-
octupole_sum(const double *oct, const vec_t *dr)
31-
{
32-
/* order in which octupoles are stored */
33-
enum { xxx = 0, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz };
34-
35-
double sum = 0.0;
36-
37-
sum += oct[xxx] * dr->x * dr->x * dr->x;
38-
sum += oct[yyy] * dr->y * dr->y * dr->y;
39-
sum += oct[zzz] * dr->z * dr->z * dr->z;
40-
sum += oct[xxy] * dr->x * dr->x * dr->y * 3.0;
41-
sum += oct[xxz] * dr->x * dr->x * dr->z * 3.0;
42-
sum += oct[xyy] * dr->x * dr->y * dr->y * 3.0;
43-
sum += oct[yyz] * dr->y * dr->y * dr->z * 3.0;
44-
sum += oct[xzz] * dr->x * dr->z * dr->z * 3.0;
45-
sum += oct[yzz] * dr->y * dr->z * dr->z * 3.0;
46-
sum += oct[xyz] * dr->x * dr->y * dr->z * 6.0;
47-
48-
return sum;
49-
}
50-
5129
static double
5230
octupole_sum_xyz(const double *oct, const vec_t *dr, size_t axis)
5331
{

‎src/pol.c‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,38 @@ get_multipole_field(const vec_t *xyz, const struct multipole_pt *mult_pt,
114114
return field;
115115
}
116116

117+
static double
118+
get_multipole_elec_potential(const vec_t *xyz, const struct multipole_pt *mult_pt,
119+
const struct swf *swf)
120+
{
121+
double elpot = 0.0;
122+
123+
vec_t dr = {
124+
xyz->x - mult_pt->x - swf->cell.x,
125+
xyz->y - mult_pt->y - swf->cell.y,
126+
xyz->z - mult_pt->z - swf->cell.z
127+
};
128+
129+
double r = vec_len(&dr);
130+
double r3 = r * r * r;
131+
double r5 = r3 * r * r;
132+
double r7 = r5 * r * r;
133+
134+
/* charge */
135+
elpot += swf->swf * mult_pt->monopole / r;
136+
137+
/* dipole */
138+
elpot += swf->swf * vec_dot(&mult_pt->dipole, &dr) / r3;
139+
140+
/* quadrupole */
141+
elpot += swf->swf * quadrupole_sum(mult_pt->quadrupole, &dr) / r5;
142+
143+
/* octupole */
144+
elpot += swf->swf * octupole_sum(mult_pt->octupole, &dr) / r7;
145+
146+
return elpot;
147+
}
148+
117149
static vec_t
118150
get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx)
119151
{
@@ -1259,3 +1291,78 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz,
12591291
*((vec_t *)field) = elec_field;
12601292
return EFP_RESULT_SUCCESS;
12611293
}
1294+
1295+
EFP_EXPORT enum efp_result
1296+
efp_get_elec_potential(struct efp *efp, size_t frag_idx, const double *xyz,
1297+
double *elec_potential)
1298+
{
1299+
assert(efp);
1300+
assert(frag_idx < efp->n_frag);
1301+
assert(xyz);
1302+
assert(elec_potential);
1303+
1304+
const struct frag *frag = efp->frags + frag_idx;
1305+
double elpot = 0.0;
1306+
1307+
for (size_t i = 0; i < efp->n_frag; i++) {
1308+
if (i == frag_idx || efp_skip_frag_pair(efp, i, frag_idx))
1309+
continue;
1310+
1311+
const struct frag *fr_i = efp->frags + i;
1312+
struct swf swf = efp_make_swf(efp, fr_i, frag, 0);
1313+
if (swf.swf == 0.0)
1314+
continue;
1315+
1316+
/* potential due to nuclei */
1317+
for (size_t j = 0; j < fr_i->n_atoms; j++) {
1318+
const struct efp_atom *at = fr_i->atoms + j;
1319+
1320+
vec_t dr = {
1321+
xyz[0] - at->x - swf.cell.x,
1322+
xyz[1] - at->y - swf.cell.y,
1323+
xyz[2] - at->z - swf.cell.z
1324+
};
1325+
1326+
double r = vec_len(&dr);
1327+
1328+
elpot += swf.swf * at->znuc / r;
1329+
}
1330+
1331+
/* potential due to multipoles */
1332+
for (size_t j = 0; j < fr_i->n_multipole_pts; j++) {
1333+
const struct multipole_pt *mpt = fr_i->multipole_pts+j;
1334+
elpot += get_multipole_elec_potential((const vec_t *)xyz, mpt, &swf);
1335+
}
1336+
1337+
/* potential due to induced dipoles */
1338+
for (size_t j = 0; j < fr_i->n_polarizable_pts; j++) {
1339+
struct polarizable_pt *pt_i = fr_i->polarizable_pts + j;
1340+
size_t idx = fr_i->polarizable_offset + j;
1341+
1342+
vec_t dr = {
1343+
xyz[0] - pt_i->x - swf.cell.x,
1344+
xyz[1] - pt_i->y - swf.cell.y,
1345+
xyz[2] - pt_i->z - swf.cell.z
1346+
};
1347+
1348+
double r = vec_len(&dr);
1349+
double r3 = r * r * r;
1350+
1351+
elpot += swf.swf * 0.5 * (vec_dot(&efp->indip[idx], &dr) + vec_dot(&efp->indipconj[idx], &dr)) / r3;
1352+
}
1353+
}
1354+
1355+
if (efp->opts.terms & EFP_TERM_AI_POL) {
1356+
/* elec potential due to nuclei from ab initio subsystem */
1357+
for (size_t i = 0; i < efp->n_ptc; i++) {
1358+
vec_t dr = vec_sub((const vec_t *)xyz, efp->ptc_xyz+i);
1359+
1360+
double r = vec_len(&dr);
1361+
1362+
elpot += efp->ptc[i] / r;
1363+
}
1364+
}
1365+
1366+
*elec_potential = elpot;
1367+
return EFP_RESULT_SUCCESS;
1368+
}

‎tests/constraint_1.in‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
run_type gtest
22
ref_energy 0.0002280010
3+
coord xyzabc
34
disp_damp tt
45
elec_damp screen
56
fraglib_path ../fraglib

‎tests/constraint_2.in‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
run_type gtest
22
ref_energy 0.0019051720
3+
coord xyzabc
34
elec_damp screen
45
disp_damp tt
56
pol_damp tt

‎tests/disp_1a.in‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
run_type gtest
22
ref_energy -0.0000989033
3+
coord xyzabc
34
terms disp
45
disp_damp tt
56
fraglib_path ../fraglib

0 commit comments

Comments
 (0)