Skip to content

Commit e52ae11

Browse files
slipchenkoslipchenko
authored
Merge pull request #12 from libefp2/full_coord
2 parents 7c2f5e0 + 32f1535 commit e52ae11

11 files changed

Lines changed: 359 additions & 29 deletions

File tree

‎efpmd/README.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,14 +59,16 @@ Default value: `sp`
5959

6060
##### Format of fragment input
6161

62-
`coord [xyzabc|points|rotmat]`
62+
`coord [xyzabc|points|rotmat|atoms]`
6363

6464
`xyzabc` - Coordinates of the center of mass and Euler angles.
6565

6666
`points` - Coordinates of three atoms for each fragment.
6767

6868
`rotmat` - Coordinates of the center of mass and rotation matrix.
6969

70+
`atoms` - Atom names and cartesian coordinates of all atoms.
71+
7072
Default value: `points`
7173

7274
See fragment input specification for more details. ATTENTION: default changed to "points" from "xyzabc" in version 1.7.3!

‎efpmd/src/common.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,9 +70,20 @@ enum ensemble_type {
7070
ENSEMBLE_TYPE_NPT
7171
};
7272

73+
/*
74+
struct atom {
75+
char *name;
76+
double x;
77+
double y;
78+
double z;
79+
};
80+
*/
81+
7382
struct frag {
7483
char *name;
7584
double coord[12];
85+
size_t n_atoms;
86+
struct efp_atom *atoms;
7687
double vel[6];
7788
bool constraint_enable;
7889
vec_t constraint_xyz;

‎efpmd/src/main.c‎

Lines changed: 33 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,10 @@
2828

2929
#include "common.h"
3030

31+
#include <stdio.h>
32+
#include <stdlib.h>
33+
#include <stdbool.h>
34+
3135
typedef void (*sim_fn_t)(struct state *);
3236

3337
void sim_sp(struct state *);
@@ -73,10 +77,12 @@ static struct cfg *make_cfg(void)
7377
cfg_add_enum(cfg, "coord", EFP_COORD_TYPE_POINTS,
7478
"xyzabc\n"
7579
"points\n"
76-
"rotmat\n",
80+
"rotmat\n"
81+
"atoms\n",
7782
(int []) { EFP_COORD_TYPE_XYZABC,
7883
EFP_COORD_TYPE_POINTS,
79-
EFP_COORD_TYPE_ROTMAT });
84+
EFP_COORD_TYPE_ROTMAT,
85+
EFP_COORD_TYPE_ATOMS});
8086

8187
cfg_add_string(cfg, "terms", "elec pol disp xr");
8288

@@ -157,7 +163,10 @@ static struct cfg *make_cfg(void)
157163
(int []) { EFP_SYMM_FRAG_FRAG,
158164
EFP_SYMM_FRAG_LIST });
159165

160-
return cfg;
166+
cfg_add_int(cfg, "update_params", 0);
167+
cfg_add_double(cfg, "update_params_cutoff", 0.0);
168+
169+
return cfg;
161170
}
162171

163172
static sim_fn_t get_sim_fn(enum run_type run_type)
@@ -275,7 +284,9 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
275284
.ligand = cfg_get_int(cfg, "ligand"),
276285
.print_pbc = cfg_get_bool(cfg, "print_pbc"),
277286
.symmetry = cfg_get_bool(cfg, "symmetry"),
278-
.symm_frag = cfg_get_enum(cfg, "symm_frag")
287+
.symm_frag = cfg_get_enum(cfg, "symm_frag"),
288+
.update_params = cfg_get_int(cfg, "update_params"),
289+
.update_params_cutoff = cfg_get_double(cfg, "update_params_cutoff")
279290
};
280291

281292
if (opts.xr_cutoff == 0.0) {
@@ -330,7 +341,16 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
330341
}
331342

332343
for (size_t i = 0; i < sys->n_frags; i++)
333-
check_fail(efp_set_frag_coordinates(efp, i, coord_type, sys->frags[i].coord));
344+
check_fail(efp_set_frag_coordinates(efp, i, coord_type, sys->frags[i].coord));
345+
346+
/*
347+
// LVS: need to coordinate this function with update_fragment() in efp.c
348+
// copying atomic coordinates of fragments
349+
// possibly add check on update_params == 1
350+
if (coord_type == EFP_COORD_TYPE_ATOMS)
351+
for (size_t i = 0; i < sys->n_frags; i++)
352+
check_fail(efp_set_frag_atoms(efp, i, sys->frags[i].n_atoms, sys->frags[i].atoms));
353+
*/
334354

335355
return (efp);
336356
}
@@ -427,7 +447,8 @@ static void convert_units(struct cfg *cfg, struct sys *sys)
427447
size_t n_convert = (size_t []) {
428448
[EFP_COORD_TYPE_XYZABC] = 3,
429449
[EFP_COORD_TYPE_POINTS] = 9,
430-
[EFP_COORD_TYPE_ROTMAT] = 3 }[cfg_get_enum(cfg, "coord")];
450+
[EFP_COORD_TYPE_ROTMAT] = 3,
451+
[EFP_COORD_TYPE_ATOMS] = 9}[cfg_get_enum(cfg, "coord")];
431452

432453
for (size_t i = 0; i < sys->n_frags; i++) {
433454
vec_scale(&sys->frags[i].constraint_xyz, 1.0 / BOHR_RADIUS);
@@ -442,8 +463,12 @@ static void convert_units(struct cfg *cfg, struct sys *sys)
442463

443464
static void sys_free(struct sys *sys)
444465
{
445-
for (size_t i = 0; i < sys->n_frags; i++)
446-
free(sys->frags[i].name);
466+
for (size_t i = 0; i < sys->n_frags; i++) {
467+
free(sys->frags[i].name);
468+
free(sys->frags[i].atoms);
469+
// for (size_t j = 0; j < sys->frags[i].n_atoms; j++)
470+
// free(sys->frags[i].atoms[j])
471+
}
447472

448473
free(sys->frags);
449474
free(sys->charges);

‎efpmd/src/parse.c‎

Lines changed: 59 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -76,20 +76,69 @@ static void parse_frag(struct stream *stream, enum efp_coord_type coord_type,
7676
int n_rows = (int []) {
7777
[EFP_COORD_TYPE_XYZABC] = 1,
7878
[EFP_COORD_TYPE_POINTS] = 3,
79-
[EFP_COORD_TYPE_ROTMAT] = 4 }[coord_type];
79+
[EFP_COORD_TYPE_ROTMAT] = 4,
80+
[EFP_COORD_TYPE_ATOMS] = 3 }[coord_type];
8081

8182
int n_cols = (int []) {
8283
[EFP_COORD_TYPE_XYZABC] = 6,
8384
[EFP_COORD_TYPE_POINTS] = 3,
84-
[EFP_COORD_TYPE_ROTMAT] = 3 }[coord_type];
85-
86-
for (int i = 0, idx = 0; i < n_rows; i++) {
87-
for (int j = 0; j < n_cols; j++, idx++)
88-
if (!efp_stream_parse_double(stream, frag->coord + idx))
89-
error("incorrect fragment coordinates format");
90-
91-
efp_stream_next_line(stream);
92-
}
85+
[EFP_COORD_TYPE_ROTMAT] = 3,
86+
[EFP_COORD_TYPE_ATOMS] = 4 }[coord_type];
87+
88+
if (coord_type == EFP_COORD_TYPE_ATOMS) {
89+
int counter = 0;
90+
while (true) {
91+
struct efp_atom atom_i;
92+
memset(&atom_i, 0, sizeof(struct efp_atom));
93+
94+
efp_stream_skip_space(stream);
95+
if (efp_strncasecmp(efp_stream_get_ptr(stream), "velocity", strlen("velocity")) == 0 ||
96+
efp_strncasecmp(efp_stream_get_ptr(stream), "constraint", strlen("constraint")) == 0 ||
97+
efp_strncasecmp(efp_stream_get_ptr(stream), "fragment", strlen("fragment")) == 0) {
98+
break;
99+
}
100+
101+
const char *ptr2 = efp_stream_get_ptr(stream);
102+
efp_stream_skip_nonspace(stream);
103+
104+
size_t len2 = efp_stream_get_ptr(stream) - ptr2;
105+
if (len2 == 0)
106+
error("incorrect fragment coordinates format: reading efp atom name");
107+
memcpy(atom_i.label, ptr2, len2);
108+
109+
efp_stream_skip_space(stream);
110+
if (!efp_stream_parse_double(stream, &atom_i.x) || !efp_stream_parse_double(stream, &atom_i.y)
111+
|| !efp_stream_parse_double(stream, &atom_i.z))
112+
error("incorrect fragment coordinates format: reading efp atom coordinates");
113+
114+
// LVS: temporary fix not to break the code; probably need to be rewritten later on
115+
if (counter < 3) {
116+
frag->coord[counter*3] = atom_i.x;
117+
frag->coord[counter*3+1] = atom_i.y;
118+
frag->coord[counter*3+2] = atom_i.z;
119+
}
120+
121+
frag->n_atoms++;
122+
frag->atoms = xrealloc(frag->atoms, frag->n_atoms * sizeof(struct efp_atom));
123+
frag->atoms[frag->n_atoms - 1] = atom_i;
124+
counter++;
125+
//printf("n_atoms = %d",frag->n_atoms);
126+
127+
efp_stream_next_line(stream);
128+
if (efp_stream_eol(stream))
129+
return;
130+
}
131+
}
132+
else {
133+
for (int i = 0, idx = 0; i < n_rows; i++) {
134+
for (int j = 0; j < n_cols; j++, idx++) {
135+
if (!efp_stream_parse_double(stream, frag->coord + idx))
136+
error("incorrect fragment coordinates format");
137+
}
138+
139+
efp_stream_next_line(stream);
140+
}
141+
}
93142

94143
efp_stream_skip_space(stream);
95144
if (efp_stream_eol(stream))

‎src/efp.c‎

Lines changed: 111 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,49 @@ set_coord_points(struct frag *frag, const double *coord)
118118
return EFP_RESULT_SUCCESS;
119119
}
120120

121+
static enum efp_result
122+
set_coord_atoms(struct frag *frag, const double *coord)
123+
{
124+
/* allow fragments with less than 3 atoms by using multipole points of
125+
* ghost atoms; multipole points have the same coordinates as atoms */
126+
if (frag->n_multipole_pts < 3) {
127+
efp_log("fragment must contain at least three atoms");
128+
return EFP_RESULT_FATAL;
129+
}
130+
131+
int natoms = frag->lib->n_atoms;
132+
133+
double ref[9] = {
134+
frag->lib->multipole_pts[0].x,
135+
frag->lib->multipole_pts[0].y,
136+
frag->lib->multipole_pts[0].z,
137+
frag->lib->multipole_pts[1].x,
138+
frag->lib->multipole_pts[1].y,
139+
frag->lib->multipole_pts[1].z,
140+
frag->lib->multipole_pts[2].x,
141+
frag->lib->multipole_pts[2].y,
142+
frag->lib->multipole_pts[2].z
143+
};
144+
145+
vec_t p1;
146+
mat_t rot1, rot2;
147+
148+
efp_points_to_matrix(coord, &rot1);
149+
efp_points_to_matrix(ref, &rot2);
150+
rot2 = mat_transpose(&rot2);
151+
frag->rotmat = mat_mat(&rot1, &rot2);
152+
p1 = mat_vec(&frag->rotmat, VEC(frag->lib->multipole_pts[0].x));
153+
154+
/* center of mass */
155+
frag->x = coord[0] - p1.x;
156+
frag->y = coord[1] - p1.y;
157+
frag->z = coord[2] - p1.z;
158+
159+
update_fragment(frag);
160+
161+
return EFP_RESULT_SUCCESS;
162+
}
163+
121164
static enum efp_result
122165
set_coord_rotmat(struct frag *frag, const double *coord)
123166
{
@@ -290,6 +333,28 @@ copy_frag(struct frag *dest, const struct frag *src)
290333
return EFP_RESULT_SUCCESS;
291334
}
292335

336+
// updates (shifts) parameters of fragment based on coordinates of fragment atoms
337+
static enum efp_result
338+
update_params(struct efp_atom *atoms, const struct frag *lib_orig, const struct frag *lib_current) {
339+
return EFP_RESULT_SUCCESS;
340+
}
341+
342+
// checks whether atoms in fragment "frag" match those in fragment "lib"
343+
static enum efp_result
344+
check_frag_atoms(struct frag *frag, const struct frag *lib) {
345+
return EFP_RESULT_SUCCESS;
346+
}
347+
348+
static enum efp_result
349+
clean_frag_atoms(struct frag *frag)
350+
{
351+
if (frag->atoms) {
352+
free(frag->atoms);
353+
frag->n_atoms = 0;
354+
}
355+
return EFP_RESULT_SUCCESS;
356+
}
357+
293358
static enum efp_result
294359
check_opts(const struct efp_opts *opts)
295360
{
@@ -885,13 +950,17 @@ efp_set_coordinates(struct efp *efp, enum efp_coord_type coord_type,
885950
case EFP_COORD_TYPE_ROTMAT:
886951
stride = 12;
887952
break;
953+
case EFP_COORD_TYPE_ATOMS:
954+
stride = 9;
955+
break;
888956
}
889957

890-
for (size_t i = 0; i < efp->n_frag; i++, coord += stride)
891-
if ((res = efp_set_frag_coordinates(efp, i, coord_type, coord))) {
958+
for (size_t i = 0; i < efp->n_frag; i++, coord += stride) {
959+
if ((res = efp_set_frag_coordinates(efp, i, coord_type, coord))) {
892960
efp_log("efp_set_coordinates() failure");
893961
return res;
894-
}
962+
}
963+
}
895964

896965
return EFP_RESULT_SUCCESS;
897966
}
@@ -915,6 +984,8 @@ efp_set_frag_coordinates(struct efp *efp, size_t frag_idx,
915984
return set_coord_points(frag, coord);
916985
case EFP_COORD_TYPE_ROTMAT:
917986
return set_coord_rotmat(frag, coord);
987+
case EFP_COORD_TYPE_ATOMS:
988+
return set_coord_atoms(frag, coord);
918989
}
919990
efp_log("efp_set_frag_coordinates() failure");
920991
assert(0);
@@ -1517,8 +1588,14 @@ efp_shutdown(struct efp *efp)
15171588
free_frag(efp->lib[i]);
15181589
free(efp->lib[i]);
15191590
}
1591+
// for the case of updated/shifted fragment parameters
1592+
for (size_t i = 0; i < efp->n_lib_current; i++) {
1593+
free_frag(efp->lib_current[i]);
1594+
free(efp->lib_current[i]);
1595+
}
15201596
free(efp->frags);
15211597
free(efp->lib);
1598+
free(efp->lib_current);
15221599
free(efp->grad);
15231600
free(efp->ptc);
15241601
free(efp->ptc_xyz);
@@ -1605,6 +1682,19 @@ efp_add_fragment(struct efp *efp, const char *name)
16051682
enum efp_result res;
16061683
struct frag *frag = efp->frags + efp->n_frag - 1;
16071684

1685+
// if update/rotate parameters
1686+
if (efp->opts.update_params == 1) {
1687+
// first, sanity check: do fragment atoms match those in the library fragment?
1688+
if (res = check_frag_atoms(frag, lib)) {
1689+
efp_log("check_frag_atoms() failure");
1690+
return res;
1691+
}
1692+
frag->rmsd = calc_rmsd(frag,lib);
1693+
if (frag->rmsd < efp->opts.update_params_cutoff) {
1694+
update_params(frag->atoms,lib, frag->lib_current);
1695+
}
1696+
}
1697+
16081698
if ((res = copy_frag(frag, lib))) {
16091699
efp_log("copy_frag() failure");
16101700
return res;
@@ -1766,6 +1856,23 @@ efp_get_frag_atoms(struct efp *efp, size_t frag_idx, size_t size,
17661856
return EFP_RESULT_SUCCESS;
17671857
}
17681858

1859+
EFP_EXPORT enum efp_result
1860+
efp_set_frag_atoms(struct efp *efp, size_t frag_idx, size_t n_atoms,
1861+
struct efp_atom *atoms)
1862+
{
1863+
struct frag *frag;
1864+
1865+
assert(efp);
1866+
assert(atoms);
1867+
assert(frag_idx < efp->n_frag);
1868+
1869+
frag = efp->frags + frag_idx;
1870+
frag->n_atoms = n_atoms;
1871+
memcpy(frag->atoms, atoms, n_atoms * sizeof(struct efp_atom));
1872+
1873+
return EFP_RESULT_SUCCESS;
1874+
}
1875+
17691876
EFP_EXPORT void
17701877
efp_torque_to_derivative(const double *euler, const double *torque,
17711878
double *deriv)
@@ -1983,3 +2090,4 @@ n_symm_frag(struct efp *efp, size_t *symm_frag) {
19832090
}
19842091
}
19852092

2093+

0 commit comments

Comments
 (0)