Skip to content

Commit 0110a85

Browse files
slipchenkoslipchenko
authored
Merge pull request #1 from libefp2/pairwise
Pairwise functionality
2 parents 15cd7ce + e76c643 commit 0110a85

65 files changed

Lines changed: 187876 additions & 6 deletions

Some content is hidden

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

‎.idea/libefp.iml‎

Lines changed: 2 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎.idea/modules.xml‎

Lines changed: 8 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎.idea/vcs.xml‎

Lines changed: 6 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎bin/cubegen.pl‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
#! /usr/bin/perl
2+
3+
# This script generates a rectangular box of fragments for EFPMD input
4+
5+
use 5.008;
6+
use strict;
7+
use warnings;
8+
9+
use Math::Trig;
10+
11+
if (scalar(@ARGV) != 6) {
12+
die <<EOF;
13+
usage: cubegen.pl <n1:n2:...> <r1:r2:...> <space> <nx> <ny> <nz>
14+
15+
<n1:n2:...> a colon separated list of fragment names
16+
<r1:r2:...> a colon separated list of fragment ratios
17+
<space> distance between the fragments
18+
<nx> number of fragments in x direction
19+
<ny> number of fragments in y direction
20+
<nz> number of fragments in z direction
21+
22+
example: cubegen.pl c2h5oh_l:h2o_l 40:60 4.0 20 20 20 > vodka.in
23+
EOF
24+
}
25+
26+
my @name = split ':', $ARGV[0];
27+
my @ratio = split ':', $ARGV[1];
28+
my $space = $ARGV[2];
29+
my $nx = $ARGV[3];
30+
my $ny = $ARGV[4];
31+
my $nz = $ARGV[5];
32+
33+
die "ratios do not match names" unless scalar(@name) == scalar(@ratio);
34+
die "positive number expected" unless $space > 0.0 && $nx > 0 && $ny > 0 && $nz > 0;
35+
36+
foreach my $i (1 .. scalar(@ratio) - 1) {
37+
$ratio[$i] += $ratio[$i - 1];
38+
}
39+
40+
foreach my $x (0 .. $nx - 1) {
41+
foreach my $y (0 .. $ny - 1) {
42+
foreach my $z (0 .. $nz - 1) {
43+
print_fragment($x, $y, $z);
44+
}
45+
}
46+
}
47+
48+
sub select_fragment {
49+
my $n = int(rand($ratio[-1]));
50+
51+
foreach my $i (0 .. scalar(@ratio) - 1) {
52+
if ($ratio[$i] > $n) {
53+
return $i;
54+
}
55+
}
56+
57+
die;
58+
}
59+
60+
sub print_fragment {
61+
my ($x, $y, $z) = @_;
62+
my $idx = select_fragment();
63+
my @xyzabc = ($space * $x,
64+
$space * $y,
65+
$space * $z,
66+
2.0 * pi * rand,
67+
pi * rand,
68+
2.0 * pi * rand);
69+
70+
print "\n", "fragment ", $name[$idx], "\n";
71+
72+
foreach (@xyzabc) {
73+
printf " %8.3f", $_;
74+
}
75+
76+
print "\n";
77+
}

‎bin/efpmd‎

392 KB
Binary file not shown.

‎bin/trajectory.pl‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
#! /usr/bin/perl
2+
3+
# This script extracts trajectory data from EFPMD output
4+
5+
use 5.008;
6+
use strict;
7+
use warnings;
8+
9+
die "usage: trajectory.pl <input>\n" if scalar(@ARGV) != 1;
10+
11+
open(FH, "<", $ARGV[0]) || die "$!";
12+
13+
while (<FH>) {
14+
next unless /GEOMETRY/;
15+
<FH>;
16+
17+
my @lines;
18+
19+
while (<FH>) {
20+
last if /^$/;
21+
push @lines, $_;
22+
}
23+
24+
print scalar(@lines), "\n";
25+
print "xyz", "\n";
26+
27+
foreach (@lines) {
28+
$_ =~ s/A[0-9]+([a-zA-Z]+)[0-9]*/$1/;
29+
print $_;
30+
}
31+
}

‎efpmd/README.md‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -215,6 +215,20 @@ Default value: `"."` (current directory)
215215
The `<path>` parameter should not contain spaces or should be in double quotes
216216
otherwise.
217217

218+
### Pairwise energy analysis
219+
220+
##### Enable/Disable pairwise analysis
221+
222+
`enable_pairwise [true|false]`
223+
224+
Default value: `false`
225+
226+
##### Specify ligand
227+
228+
`ligand [integer]`
229+
230+
Default value: `0` (the first fragment in the system)
231+
218232
### Periodic Boundary Conditions (PBC)
219233

220234
##### Enable/Disable PBC

‎efpmd/src/common.c‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -168,6 +168,7 @@ void print_energy(struct state *state)
168168
msg("\n\n");
169169
}
170170

171+
171172
void print_gradient(struct state *state)
172173
{
173174
size_t n_frags;
@@ -270,6 +271,97 @@ void print_matrix(size_t rows, size_t cols, const double *mat)
270271
}
271272
}
272273

274+
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;
301+
}
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;
348+
}
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+
}
362+
}
363+
364+
273365
vec_t box_from_str(const char *str)
274366
{
275367
vec_t box;
@@ -280,3 +372,4 @@ vec_t box_from_str(const char *str)
280372
vec_scale(&box, 1.0 / BOHR_RADIUS);
281373
return box;
282374
}
375+

‎efpmd/src/common.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,8 +59,9 @@ enum run_type {
5959
RUN_TYPE_OPT,
6060
RUN_TYPE_MD,
6161
RUN_TYPE_EFIELD,
62-
RUN_TYPE_GTEST
62+
RUN_TYPE_GTEST,
6363
};
64+
//RUN_TYPE_PAIRWISE
6465

6566
enum ensemble_type {
6667
ENSEMBLE_TYPE_NVE,
@@ -113,6 +114,7 @@ void print_fragment(const char *, const double *, const double *);
113114
void print_charge(double, double, double, double);
114115
void print_vector(size_t, const double *);
115116
void print_matrix(size_t, size_t, const double *);
117+
void print_pair_energy(struct state *);
116118

117119
void check_fail(enum efp_result);
118120
void compute_energy(struct state *, bool);
@@ -121,4 +123,5 @@ vec_t box_from_str(const char *);
121123
int efp_strcasecmp(const char *, const char *);
122124
int efp_strncasecmp(const char *, const char *, size_t);
123125

126+
124127
#endif /* EFPMD_COMMON_H */

‎efpmd/src/main.c‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ 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 *);
4041

4142
#define USAGE_STRING \
4243
"usage: efpmd [-d | -v | -h | input]\n" \
@@ -62,7 +63,9 @@ static struct cfg *make_cfg(void)
6263
RUN_TYPE_OPT,
6364
RUN_TYPE_MD,
6465
RUN_TYPE_EFIELD,
65-
RUN_TYPE_GTEST });
66+
RUN_TYPE_GTEST});
67+
/* "pairwise\n",
68+
RUN_TYPE_PAIRWISE */
6669

6770
cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_XYZABC,
6871
"xyzabc\n"
@@ -139,6 +142,10 @@ static struct cfg *make_cfg(void)
139142
cfg_add_double(cfg, "thermostat_tau", 1.0e3);
140143
cfg_add_double(cfg, "barostat_tau", 1.0e4);
141144

145+
cfg_add_int(cfg, "ligand", 0);
146+
cfg_add_bool(cfg, "enable_pairwise", false);
147+
148+
142149
return cfg;
143150
}
144151

@@ -160,6 +167,8 @@ static sim_fn_t get_sim_fn(enum run_type run_type)
160167
case RUN_TYPE_GTEST:
161168
return sim_gtest;
162169
}
170+
// case RUN_TYPE_PAIRWISE:
171+
// return sim_ pairwise;
163172
assert(0);
164173
}
165174

@@ -247,7 +256,9 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
247256
.pol_driver = cfg_get_enum(cfg, "pol_driver"),
248257
.enable_pbc = cfg_get_bool(cfg, "enable_pbc"),
249258
.enable_cutoff = cfg_get_bool(cfg, "enable_cutoff"),
250-
.swf_cutoff = cfg_get_double(cfg, "swf_cutoff")
259+
.swf_cutoff = cfg_get_double(cfg, "swf_cutoff"),
260+
.enable_pairwise = cfg_get_bool(cfg, "enable_pairwise"),
261+
.ligand = cfg_get_int(cfg, "ligand")
251262
};
252263

253264
enum efp_coord_type coord_type = cfg_get_enum(cfg, "coord");
@@ -310,6 +321,7 @@ static void state_init(struct state *state, const struct cfg *cfg, const struct
310321
state->grad = xcalloc(sys->n_frags * 6 + sys->n_charges * 3, sizeof(double));
311322
state->ff = NULL;
312323

324+
313325
if (cfg_get_bool(cfg, "enable_ff")) {
314326
if ((state->ff = ff_create()) == NULL)
315327
error("cannot create ff object");

0 commit comments

Comments
 (0)