1 /*-
2 * Copyright (c) 2012-2017 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 <stdlib.h>
28
29 #include "clapack.h"
30 #include "private.h"
31
32 double efp_get_pol_damp_tt(double, double, double);
33 enum efp_result efp_compute_id_direct(struct efp *);
34
35 static void
copy_matrix(double * dst,size_t n,size_t off_i,size_t off_j,const mat_t * m)36 copy_matrix(double *dst, size_t n, size_t off_i, size_t off_j, const mat_t *m)
37 {
38 dst[n * (3 * off_i + 0) + 3 * off_j + 0] = m->xx;
39 dst[n * (3 * off_i + 0) + 3 * off_j + 1] = m->xy;
40 dst[n * (3 * off_i + 0) + 3 * off_j + 2] = m->xz;
41 dst[n * (3 * off_i + 1) + 3 * off_j + 0] = m->yx;
42 dst[n * (3 * off_i + 1) + 3 * off_j + 1] = m->yy;
43 dst[n * (3 * off_i + 1) + 3 * off_j + 2] = m->yz;
44 dst[n * (3 * off_i + 2) + 3 * off_j + 0] = m->zx;
45 dst[n * (3 * off_i + 2) + 3 * off_j + 1] = m->zy;
46 dst[n * (3 * off_i + 2) + 3 * off_j + 2] = m->zz;
47 }
48
49 static void
transpose_matrix(double * m,size_t n)50 transpose_matrix(double *m, size_t n)
51 {
52 for (size_t i = 0; i < n; i++)
53 for (size_t j = i + 1; j < n; j++) {
54 double t = m[n * i + j];
55 m[n * i + j] = m[n * j + i];
56 m[n * j + i] = t;
57 }
58 }
59
60 static mat_t
get_int_mat(const struct efp * efp,size_t i,size_t j,size_t ii,size_t jj)61 get_int_mat(const struct efp *efp, size_t i, size_t j, size_t ii, size_t jj)
62 {
63 mat_t m;
64 const struct frag *fr_i = efp->frags + i;
65 const struct frag *fr_j = efp->frags + j;
66 const struct polarizable_pt *pt_i = fr_i->polarizable_pts + ii;
67 const struct polarizable_pt *pt_j = fr_j->polarizable_pts + jj;
68 struct swf swf = efp_make_swf(efp, fr_i, fr_j);
69
70 vec_t dr = {
71 pt_j->x - pt_i->x - swf.cell.x,
72 pt_j->y - pt_i->y - swf.cell.y,
73 pt_j->z - pt_i->z - swf.cell.z
74 };
75
76 double p1 = 1.0;
77 double r = vec_len(&dr);
78 double r3 = r * r * r;
79 double r5 = r3 * r * r;
80
81 if (efp->opts.pol_damp == EFP_POL_DAMP_TT)
82 p1 = efp_get_pol_damp_tt(r, fr_i->pol_damp, fr_j->pol_damp);
83
84 m.xx = swf.swf * p1 * (3.0 * dr.x * dr.x / r5 - 1.0 / r3);
85 m.xy = swf.swf * p1 * 3.0 * dr.x * dr.y / r5;
86 m.xz = swf.swf * p1 * 3.0 * dr.x * dr.z / r5;
87 m.yx = swf.swf * p1 * 3.0 * dr.y * dr.x / r5;
88 m.yy = swf.swf * p1 * (3.0 * dr.y * dr.y / r5 - 1.0 / r3);
89 m.yz = swf.swf * p1 * 3.0 * dr.y * dr.z / r5;
90 m.zx = swf.swf * p1 * 3.0 * dr.z * dr.x / r5;
91 m.zy = swf.swf * p1 * 3.0 * dr.z * dr.y / r5;
92 m.zz = swf.swf * p1 * (3.0 * dr.z * dr.z / r5 - 1.0 / r3);
93
94 return m;
95 }
96
97 static void
compute_lhs(const struct efp * efp,double * c,int conj)98 compute_lhs(const struct efp *efp, double *c, int conj)
99 {
100 size_t i, ii, j, jj, offset_i, offset_j;
101 size_t n = 3 * efp->n_polarizable_pts;
102
103 for (i = 0, offset_i = 0; i < efp->n_frag; i++) {
104 for (ii = 0; ii < efp->frags[i].n_polarizable_pts; ii++, offset_i++) {
105 for (j = 0, offset_j = 0; j < efp->n_frag; j++) {
106 for (jj = 0; jj < efp->frags[j].n_polarizable_pts; jj++, offset_j++) {
107 if (i == j) {
108 if (ii == jj) {
109 copy_matrix(c, n, offset_i, offset_j,
110 &mat_identity);
111 } else {
112 copy_matrix(c, n, offset_i, offset_j,
113 &mat_zero);
114 }
115 continue;
116 }
117 const struct polarizable_pt *pt_i =
118 efp->frags[i].polarizable_pts + ii;
119 mat_t m = get_int_mat(efp, i, j, ii, jj);
120
121 if (conj)
122 m = mat_trans_mat(&pt_i->tensor, &m);
123 else
124 m = mat_mat(&pt_i->tensor, &m);
125
126 mat_negate(&m);
127 copy_matrix(c, n, offset_i, offset_j, &m);
128 }}}}
129 }
130
131 static void
compute_rhs(const struct efp * efp,vec_t * id,int conj)132 compute_rhs(const struct efp *efp, vec_t *id, int conj)
133 {
134 for (size_t i = 0, idx = 0; i < efp->n_frag; i++) {
135 const struct frag *frag = efp->frags + i;
136
137 for (size_t j = 0; j < frag->n_polarizable_pts; j++, idx++) {
138 const struct polarizable_pt *pt =
139 frag->polarizable_pts + j;
140 vec_t field = vec_add(&pt->elec_field,
141 &pt->elec_field_wf);
142
143 if (conj)
144 id[idx] = mat_trans_vec(&pt->tensor, &field);
145 else
146 id[idx] = mat_vec(&pt->tensor, &field);
147 }
148 }
149 }
150
151 enum efp_result
efp_compute_id_direct(struct efp * efp)152 efp_compute_id_direct(struct efp *efp)
153 {
154 double *c;
155 size_t n;
156 fortranint_t *ipiv;
157 enum efp_result res;
158
159 n = 3 * efp->n_polarizable_pts;
160 c = (double *)calloc(n * n, sizeof *c);
161 ipiv = (fortranint_t *)calloc(n, sizeof *ipiv);
162
163 if (c == NULL || ipiv == NULL) {
164 res = EFP_RESULT_NO_MEMORY;
165 goto error;
166 }
167
168 /* induced dipoles */
169 compute_lhs(efp, c, 0);
170 compute_rhs(efp, efp->indip, 0);
171 transpose_matrix(c, n);
172
173 if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv,
174 (double *)efp->indip, (fortranint_t)n) != 0) {
175 efp_log("dgesv: error solving for induced dipoles");
176 res = EFP_RESULT_FATAL;
177 goto error;
178 }
179
180 /* conjugate induced dipoles */
181 compute_lhs(efp, c, 1);
182 compute_rhs(efp, efp->indipconj, 1);
183 transpose_matrix(c, n);
184
185 if (efp_dgesv((fortranint_t)n, 1, c, (fortranint_t)n, ipiv,
186 (double *)efp->indipconj, (fortranint_t)n) != 0) {
187 efp_log("dgesv: error solving for conjugate induced dipoles");
188 res = EFP_RESULT_FATAL;
189 goto error;
190 }
191 res = EFP_RESULT_SUCCESS;
192 error:
193 free(c);
194 free(ipiv);
195 return res;
196 }
197