1 /*-
2 * Copyright (c) 2012-2015 Ilya Kaliman
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27 #include "common.h"
28
29 /* current coordinates from efp struct are used */
compute_energy(struct state * state,bool do_grad)30 void compute_energy(struct state *state, bool do_grad)
31 {
32 struct efp_atom *atoms;
33 struct efp_energy efp_energy;
34 double xyz[3], xyzabc[6], *grad;
35 size_t ifrag, nfrag, iatom, natom;
36 int itotal;
37
38 /* EFP part */
39 check_fail(efp_compute(state->efp, do_grad));
40 check_fail(efp_get_energy(state->efp, &efp_energy));
41 check_fail(efp_get_frag_count(state->efp, &nfrag));
42
43 if (do_grad) {
44 check_fail(efp_get_gradient(state->efp, state->grad));
45 check_fail(efp_get_point_charge_gradient(state->efp,
46 state->grad + 6 * nfrag));
47 }
48
49 state->energy = efp_energy.total;
50
51 /* constraints */
52 for (ifrag = 0; ifrag < nfrag; ifrag++) {
53 const struct frag *frag = state->sys->frags + ifrag;
54
55 check_fail(efp_get_frag_xyzabc(state->efp, ifrag, xyzabc));
56
57 if (frag->constraint_enable) {
58 double dr2, drx, dry, drz;
59
60 drx = xyzabc[0] - frag->constraint_xyz.x;
61 dry = xyzabc[1] - frag->constraint_xyz.y;
62 drz = xyzabc[2] - frag->constraint_xyz.z;
63
64 dr2 = drx * drx + dry * dry + drz * drz;
65 state->energy += 0.5 * frag->constraint_k * dr2;
66
67 if (do_grad) {
68 grad = state->grad + 6 * ifrag;
69 grad[0] += frag->constraint_k * drx;
70 grad[1] += frag->constraint_k * dry;
71 grad[2] += frag->constraint_k * drz;
72 }
73 }
74 }
75
76 /* MM force field part */
77 if (state->ff == NULL)
78 return;
79
80 for (ifrag = 0, itotal = 0; ifrag < nfrag; ifrag++) {
81 check_fail(efp_get_frag_atom_count(state->efp, ifrag, &natom));
82 atoms = xmalloc(natom * sizeof(struct efp_atom));
83 check_fail(efp_get_frag_atoms(state->efp, ifrag, natom, atoms));
84
85 for (iatom = 0; iatom < natom; iatom++, itotal++)
86 ff_set_atom_xyz(state->ff, itotal, &atoms[iatom].x);
87
88 free(atoms);
89 }
90
91 ff_compute(state->ff, do_grad);
92
93 if (do_grad) {
94 for (ifrag = 0, itotal = 0, grad = state->grad; ifrag < nfrag; ifrag++, grad += 6) {
95 check_fail(efp_get_frag_xyzabc(state->efp, ifrag, xyzabc));
96 check_fail(efp_get_frag_atom_count(state->efp, ifrag, &natom));
97 atoms = xmalloc(natom * sizeof(struct efp_atom));
98 check_fail(efp_get_frag_atoms(state->efp, ifrag, natom, atoms));
99
100 for (iatom = 0; iatom < natom; iatom++, itotal++) {
101 ff_get_atom_gradient(state->ff, itotal, xyz);
102
103 grad[0] += xyz[0];
104 grad[1] += xyz[1];
105 grad[2] += xyz[2];
106
107 grad[3] += (atoms[iatom].y - xyzabc[1]) * xyz[2] -
108 (atoms[iatom].z - xyzabc[2]) * xyz[1];
109 grad[4] += (atoms[iatom].z - xyzabc[2]) * xyz[0] -
110 (atoms[iatom].x - xyzabc[0]) * xyz[2];
111 grad[5] += (atoms[iatom].x - xyzabc[0]) * xyz[1] -
112 (atoms[iatom].y - xyzabc[1]) * xyz[0];
113 }
114
115 free(atoms);
116 }
117 }
118
119 state->energy += ff_get_energy(state->ff);
120 }
121