Skip to content

Commit 01145a7

Browse files
slipchenkoslipchenko
committed
nonorthogonal PBC and PBC print-out
1 parent 3f4ff27 commit 01145a7

4 files changed

Lines changed: 146 additions & 96 deletions

File tree

‎efpmd/src/common.c‎

Lines changed: 137 additions & 87 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,61 @@ void print_geometry(struct efp *efp)
140140
msg("\n\n");
141141
}
142142

143+
void print_geometry_pbc(struct efp *efp, int ligand)
144+
{
145+
size_t n_frags;
146+
check_fail(efp_get_frag_count(efp, &n_frags));
147+
148+
msg(" GEOMETRY IN PBC CELL (ANGSTROMS)\n\n");
149+
150+
double bx[6];
151+
152+
check_fail(efp_get_periodic_box(efp, bx));
153+
six_t box = {bx[0],bx[1],bx[2],bx[3],bx[4],bx[5]};
154+
//six_t box = {10.66, 12.03, 10.872, 90.0, 115.83, 90.0};
155+
156+
double lig_com[6];
157+
check_fail(efp_get_frag_xyzabc(efp,ligand,lig_com));
158+
for (size_t i = 0; i < n_frags; i++) {
159+
160+
vec_t cell = {0.0,0.0,0.0};
161+
if (i != ligand) {
162+
double frag_com[6];
163+
check_fail(efp_get_frag_xyzabc(efp, i, frag_com));
164+
vec_t dr = {frag_com[0] - lig_com[0], frag_com[1] - lig_com[1], frag_com[2] - lig_com[2]};
165+
166+
if (box.a == 90.0 && box.b == 90.0 && box.c == 90.0) {
167+
cell.x = box.x * round(dr.x / box.x);
168+
cell.y = box.y * round(dr.y / box.y);
169+
cell.z = box.z * round(dr.z / box.z);
170+
//dr = vec_sub(&dr, &cell);
171+
} else {
172+
cart_to_frac(box, &dr);
173+
cell.x = round(dr.x);
174+
cell.y = round(dr.y);
175+
cell.z = round(dr.z);
176+
//dr = vec_sub(&dr, &cell);
177+
frac_to_cart(box, &cell);
178+
}
179+
}
180+
size_t n_atoms;
181+
check_fail(efp_get_frag_atom_count(efp, i, &n_atoms));
182+
183+
struct efp_atom atoms[n_atoms];
184+
check_fail(efp_get_frag_atoms(efp, i, n_atoms, atoms));
185+
186+
for (size_t a = 0; a < n_atoms; a++) {
187+
double x = (atoms[a].x - cell.x)* BOHR_RADIUS;
188+
double y = (atoms[a].y - cell.y)* BOHR_RADIUS;
189+
double z = (atoms[a].z - cell.z)* BOHR_RADIUS;
190+
191+
msg("%-16s %12.6lf %12.6lf %12.6lf\n", atoms[a].label, x, y, z);
192+
}
193+
}
194+
195+
msg("\n\n");
196+
}
197+
143198
void print_energy(struct state *state)
144199
{
145200
struct efp_energy energy;
@@ -272,97 +327,92 @@ void print_matrix(size_t rows, size_t cols, const double *mat)
272327
}
273328

274329
void print_pair_energy(struct state *state){
275-
276-
if (cfg_get_bool(state->cfg, "enable_pairwise")) {
277-
msg(" ------ PAIRWISE ENERGY ANALYSIS FOLLOWS ------------------\n\n");
278-
size_t n_frags;
279-
check_fail(efp_get_frag_count(state->efp, &n_frags));
280-
double coord[6 * n_frags];
281-
check_fail(efp_get_coordinates(state->efp, coord));
282-
283-
struct efp_energy *energies;
284-
energies = xmalloc(n_frags * sizeof(struct efp_energy));
285-
check_fail(efp_get_pairwise_energy(state->efp, energies));
286-
287-
char ligand[32];
288-
size_t lig_atoms;
289-
290-
size_t ligand_index = cfg_get_int(state->cfg, "ligand");
291-
check_fail(efp_get_frag_name(state->efp, ligand_index, sizeof(ligand),ligand));
292-
check_fail(efp_get_frag_atom_count(state->efp, ligand_index, &lig_atoms));
293-
struct efp_atom latoms[lig_atoms];
294-
check_fail(efp_get_frag_atoms(state->efp, ligand_index, lig_atoms, latoms));
295-
296-
char frag_name[32];
297-
size_t frag_atoms;
298-
double lattice_energy[6];
299-
for (size_t j=0; j<6; j++){
300-
lattice_energy[j]=0.0;
330+
msg(" ------ PAIRWISE ENERGY ANALYSIS FOLLOWS ------------------\n\n");
331+
size_t n_frags;
332+
check_fail(efp_get_frag_count(state->efp, &n_frags));
333+
double coord[6 * n_frags];
334+
check_fail(efp_get_coordinates(state->efp, coord));
335+
336+
struct efp_energy *energies;
337+
energies = xmalloc(n_frags * sizeof(struct efp_energy));
338+
check_fail(efp_get_pairwise_energy(state->efp, energies));
339+
340+
char ligand[32];
341+
size_t lig_atoms;
342+
343+
size_t ligand_index = cfg_get_int(state->cfg, "ligand");
344+
check_fail(efp_get_frag_name(state->efp, ligand_index, sizeof(ligand),ligand));
345+
check_fail(efp_get_frag_atom_count(state->efp, ligand_index, &lig_atoms));
346+
struct efp_atom latoms[lig_atoms];
347+
check_fail(efp_get_frag_atoms(state->efp, ligand_index, lig_atoms, latoms));
348+
349+
char frag_name[32];
350+
size_t frag_atoms;
351+
double lattice_energy[6];
352+
for (size_t j=0; j<6; j++){
353+
lattice_energy[j]=0.0;
354+
}
355+
for (size_t i=0; i <n_frags; i++){
356+
check_fail(efp_get_frag_name(state->efp, i, sizeof(frag_name),frag_name));
357+
check_fail(efp_get_frag_atom_count(state->efp, i, &frag_atoms));
358+
359+
struct efp_atom atoms[frag_atoms];
360+
check_fail(efp_get_frag_atoms(state->efp, i, frag_atoms, atoms));
361+
362+
msg(" PAIRWISE ENERGY BETWEEN FRAGMENT %zu (%s) AND LIGAND %zu (%s) \n", i, frag_name, ligand_index, ligand);
363+
msg("fragment %s\n", frag_name);
364+
for (size_t a = 0; a < frag_atoms; a++) {
365+
double x = atoms[a].x * BOHR_RADIUS;
366+
double y = atoms[a].y * BOHR_RADIUS;
367+
double z = atoms[a].z * BOHR_RADIUS;
368+
msg(" %-16s %12.6lf %12.6lf %12.6lf\n", atoms[a].label, x, y, z);
301369
}
302-
for (size_t i=0; i <n_frags; i++){
303-
check_fail(efp_get_frag_name(state->efp, i, sizeof(frag_name),frag_name));
304-
check_fail(efp_get_frag_atom_count(state->efp, i, &frag_atoms));
305-
306-
struct efp_atom atoms[frag_atoms];
307-
check_fail(efp_get_frag_atoms(state->efp, i, frag_atoms, atoms));
308-
309-
msg(" PAIRWISE ENERGY BETWEEN FRAGMENT %zu (%s) AND LIGAND %zu (%s) \n", i, frag_name, ligand_index, ligand);
310-
msg("fragment %s\n", frag_name);
311-
for (size_t a = 0; a < frag_atoms; a++) {
312-
double x = atoms[a].x * BOHR_RADIUS;
313-
double y = atoms[a].y * BOHR_RADIUS;
314-
double z = atoms[a].z * BOHR_RADIUS;
315-
msg(" %-16s %12.6lf %12.6lf %12.6lf\n", atoms[a].label, x, y, z);
316-
}
317-
msg("\n");
318-
319-
msg("fragment %s\n", ligand);
320-
for (size_t a = 0; a < lig_atoms; a++) {
321-
double x = latoms[a].x * BOHR_RADIUS;
322-
double y = latoms[a].y * BOHR_RADIUS;
323-
double z = latoms[a].z * BOHR_RADIUS;
324-
msg(" %-16s %12.6lf %12.6lf %12.6lf\n", latoms[a].label, x, y, z);
325-
}
326-
msg("\n");
327-
328-
msg(" PAIRWISE ENERGY COMPONENTS (ATOMIC UNITS)\n");
329-
msg("%40s %16.10lf\n", "PAIRWISE ELECTROSTATIC ENERGY", energies[i].electrostatic);
330-
msg("%40s %16.10lf\n", "PAIRWISE POLARIZATION ENERGY", energies[i].polarization);
331-
msg("%40s %16.10lf\n", "PAIRWISE DISPERSION ENERGY", energies[i].dispersion);
332-
msg("%40s %16.10lf\n", "PAIRWISE EXCHANGE REPULSION ENERGY",
333-
energies[i].exchange_repulsion);
334-
msg("%40s %16.10lf\n", "PAIRWISE CHARGE PENETRATION ENERGY",
335-
energies[i].charge_penetration);
336-
337-
energies[i].total = energies[i].electrostatic + energies[i].polarization + energies[i].dispersion
338-
+ energies[i].exchange_repulsion + energies[i].charge_penetration;
339-
msg("%40s %16.10lf\n", "PAIRWISE TOTAL ENERGY", energies[i].total);
340-
msg("\n ---------------------------------------------------------\n");
341-
342-
lattice_energy[0] = lattice_energy[0] + energies[i].electrostatic;
343-
lattice_energy[1] = lattice_energy[1] + energies[i].polarization;
344-
lattice_energy[2] = lattice_energy[2] + energies[i].dispersion;
345-
lattice_energy[3] = lattice_energy[3] + energies[i].exchange_repulsion;
346-
lattice_energy[4] = lattice_energy[4] + energies[i].charge_penetration;
347-
lattice_energy[5] = lattice_energy[5] + energies[i].total;
370+
msg("\n");
371+
372+
msg("fragment %s\n", ligand);
373+
for (size_t a = 0; a < lig_atoms; a++) {
374+
double x = latoms[a].x * BOHR_RADIUS;
375+
double y = latoms[a].y * BOHR_RADIUS;
376+
double z = latoms[a].z * BOHR_RADIUS;
377+
msg(" %-16s %12.6lf %12.6lf %12.6lf\n", latoms[a].label, x, y, z);
348378
}
349-
free(energies);
350-
351-
msg(" LATTICE ENERGY COMPONENTS (ATOMIC UNITS)\n");
352-
msg("%40s %16.10lf\n", "LATTICE ELECTROSTATIC ENERGY", lattice_energy[0]);
353-
msg("%40s %16.10lf\n", "LATTICE POLARIZATION ENERGY", lattice_energy[1]);
354-
msg("%40s %16.10lf\n", "LATTICE DISPERSION ENERGY", lattice_energy[2]);
355-
msg("%40s %16.10lf\n", "LATTICE EXCHANGE REPULSION ENERGY", lattice_energy[3]);
356-
msg("%40s %16.10lf\n", "LATTICE CHARGE PENETRATION ENERGY", lattice_energy[4]);
357-
msg("%40s %16.10lf\n", "LATTICE TOTAL ENERGY", lattice_energy[5]);
358-
msg("\n\n");
359-
360-
msg(" ------ PAIRWISE ENERGY ANALYSIS COMPLETED --------------\n\n");
361-
}
379+
msg("\n");
380+
381+
msg(" PAIRWISE ENERGY COMPONENTS (ATOMIC UNITS)\n");
382+
msg("%40s %16.10lf\n", "PAIRWISE ELECTROSTATIC ENERGY", energies[i].electrostatic);
383+
msg("%40s %16.10lf\n", "PAIRWISE POLARIZATION ENERGY", energies[i].polarization);
384+
msg("%40s %16.10lf\n", "PAIRWISE DISPERSION ENERGY", energies[i].dispersion);
385+
msg("%40s %16.10lf\n", "PAIRWISE EXCHANGE REPULSION ENERGY",
386+
energies[i].exchange_repulsion);
387+
msg("%40s %16.10lf\n", "PAIRWISE CHARGE PENETRATION ENERGY",
388+
energies[i].charge_penetration);
389+
390+
energies[i].total = energies[i].electrostatic + energies[i].polarization + energies[i].dispersion
391+
+ energies[i].exchange_repulsion + energies[i].charge_penetration;
392+
msg("%40s %16.10lf\n", "PAIRWISE TOTAL ENERGY", energies[i].total);
393+
msg("\n ---------------------------------------------------------\n");
394+
395+
lattice_energy[0] = lattice_energy[0] + energies[i].electrostatic;
396+
lattice_energy[1] = lattice_energy[1] + energies[i].polarization;
397+
lattice_energy[2] = lattice_energy[2] + energies[i].dispersion;
398+
lattice_energy[3] = lattice_energy[3] + energies[i].exchange_repulsion;
399+
lattice_energy[4] = lattice_energy[4] + energies[i].charge_penetration;
400+
lattice_energy[5] = lattice_energy[5] + energies[i].total;
401+
}
402+
free(energies);
403+
404+
msg(" LATTICE ENERGY COMPONENTS (ATOMIC UNITS)\n");
405+
msg("%40s %16.10lf\n", "LATTICE ELECTROSTATIC ENERGY", lattice_energy[0]);
406+
msg("%40s %16.10lf\n", "LATTICE POLARIZATION ENERGY", lattice_energy[1]);
407+
msg("%40s %16.10lf\n", "LATTICE DISPERSION ENERGY", lattice_energy[2]);
408+
msg("%40s %16.10lf\n", "LATTICE EXCHANGE REPULSION ENERGY", lattice_energy[3]);
409+
msg("%40s %16.10lf\n", "LATTICE CHARGE PENETRATION ENERGY", lattice_energy[4]);
410+
msg("%40s %16.10lf\n", "LATTICE TOTAL ENERGY", lattice_energy[5]);
411+
msg("\n\n");
412+
413+
msg(" ------ PAIRWISE ENERGY ANALYSIS COMPLETED --------------\n\n");
362414
}
363415

364-
365-
366416
six_t box_from_str(const char *str)
367417
{
368418
six_t box;

‎efpmd/src/common.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,6 @@ enum run_type {
6161
RUN_TYPE_EFIELD,
6262
RUN_TYPE_GTEST,
6363
};
64-
//RUN_TYPE_PAIRWISE
6564

6665
enum ensemble_type {
6766
ENSEMBLE_TYPE_NVE,
@@ -108,6 +107,7 @@ void *xrealloc(void *, size_t);
108107

109108
void print_vec(const double *);
110109
void print_geometry(struct efp *);
110+
void print_geometry_pbc(struct efp *, int);
111111
void print_energy(struct state *);
112112
void print_gradient(struct state *);
113113
void print_fragment(const char *, const double *, const double *);

‎efpmd/src/main.c‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,6 @@ void sim_opt(struct state *);
3737
void sim_md(struct state *);
3838
void sim_efield(struct state *);
3939
void sim_gtest(struct state *);
40-
// void sim_pairwise(struct state *);
4140

4241
#define USAGE_STRING \
4342
"usage: efpmd [-d | -v | -h | input]\n" \
@@ -64,8 +63,6 @@ static struct cfg *make_cfg(void)
6463
RUN_TYPE_MD,
6564
RUN_TYPE_EFIELD,
6665
RUN_TYPE_GTEST});
67-
/* "pairwise\n",
68-
RUN_TYPE_PAIRWISE */
6966

7067
cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
7168
"xyzabc\n"
@@ -143,7 +140,8 @@ static struct cfg *make_cfg(void)
143140
cfg_add_double(cfg, "barostat_tau", 1.0e4);
144141

145142
cfg_add_int(cfg, "ligand", 0);
146-
cfg_add_bool(cfg, "enable_pairwise", false);
143+
cfg_add_bool(cfg, "enable_pairwise", false);
144+
cfg_add_bool(cfg, "print_pbc", false);
147145

148146

149147
return cfg;
@@ -167,8 +165,6 @@ static sim_fn_t get_sim_fn(enum run_type run_type)
167165
case RUN_TYPE_GTEST:
168166
return sim_gtest;
169167
}
170-
// case RUN_TYPE_PAIRWISE:
171-
// return sim_ pairwise;
172168
assert(0);
173169
}
174170

@@ -258,7 +254,8 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
258254
.enable_cutoff = cfg_get_bool(cfg, "enable_cutoff"),
259255
.swf_cutoff = cfg_get_double(cfg, "swf_cutoff"),
260256
.enable_pairwise = cfg_get_bool(cfg, "enable_pairwise"),
261-
.ligand = cfg_get_int(cfg, "ligand")
257+
.ligand = cfg_get_int(cfg, "ligand"),
258+
.print_pbc = cfg_get_bool(cfg, "print_pbc")
262259
};
263260

264261
enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord");

‎efpmd/src/sp.c‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,11 @@ void sim_sp(struct state *state)
3333
msg("SINGLE POINT ENERGY JOB\n\n\n");
3434

3535
print_geometry(state->efp);
36+
if (cfg_get_bool(state->cfg, "print_pbc") && cfg_get_bool(state->cfg, "enable_pbc"))
37+
print_geometry_pbc(state->efp,cfg_get_int(state->cfg, "ligand"));
3638
compute_energy(state, false);
37-
print_pair_energy(state);
39+
if (cfg_get_bool(state->cfg, "enable_pairwise"))
40+
print_pair_energy(state);
3841
print_energy(state);
3942

4043
msg("SINGLE POINT ENERGY JOB COMPLETED SUCCESSFULLY\n");

0 commit comments

Comments
 (0)