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 "private.h"
28
29 static const double weights[] = {
30 0.72086099022968040154E-02, 0.17697067815034886394E-01,
31 0.30660908596251749739E-01, 0.48381293256249884995E-01,
32 0.74878830420650517080E-01, 0.11806515901361630228E+00,
33 0.19535413832209084204E+00, 0.35055692324483221824E+00,
34 0.71577113554429568336E+00, 0.18140975997632396972E+01,
35 0.69792344511487082324E+01, 0.83248093882965845391E+02
36 };
37
38 static double
get_damp_tt(double r)39 get_damp_tt(double r)
40 {
41 static const double a = 1.5; /* Tang-Toennies damping parameter */
42
43 double ra = r * a;
44 double ra2 = ra * ra;
45 double ra3 = ra2 * ra;
46 double ra4 = ra3 * ra;
47 double ra5 = ra4 * ra;
48 double ra6 = ra5 * ra;
49
50 return 1.0 - exp(-ra) * (1.0 + ra + ra2 / 2.0 + ra3 / 6.0 +
51 ra4 / 24.0 + ra5 / 120.0 + ra6 / 720.0);
52 }
53
54 static double
get_damp_tt_grad(double r)55 get_damp_tt_grad(double r)
56 {
57 static const double a = 1.5; /* Tang-Toennies damping parameter */
58
59 double ra = r * a;
60 double ra2 = ra * ra;
61 double ra6 = ra2 * ra2 * ra2;
62
63 return a * exp(-ra) * ra6 / 720.0;
64 }
65
66 static double
disp_tt(struct efp * efp,size_t fr_i_idx,size_t fr_j_idx,size_t pt_i_idx,size_t pt_j_idx,double sum,const struct swf * swf)67 disp_tt(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, size_t pt_i_idx,
68 size_t pt_j_idx, double sum, const struct swf *swf)
69 {
70 const struct frag *fr_i = efp->frags + fr_i_idx;
71 const struct frag *fr_j = efp->frags + fr_j_idx;
72
73 const struct dynamic_polarizable_pt *pt_i =
74 fr_i->dynamic_polarizable_pts + pt_i_idx;
75 const struct dynamic_polarizable_pt *pt_j =
76 fr_j->dynamic_polarizable_pts + pt_j_idx;
77
78 vec_t dr = {
79 pt_j->x - pt_i->x - swf->cell.x,
80 pt_j->y - pt_i->y - swf->cell.y,
81 pt_j->z - pt_i->z - swf->cell.z
82 };
83
84 double r = vec_len(&dr);
85 double r2 = r * r;
86 double r6 = r2 * r2 * r2;
87
88 double damp = get_damp_tt(r);
89 double energy = -4.0 / 3.0 * sum * damp / r6;
90
91 if (efp->do_gradient) {
92 double gdamp = get_damp_tt_grad(r);
93 double g = 4.0 / 3.0 * sum * (gdamp / r - 6.0 * damp / r2) / r6;
94
95 vec_t force = {
96 g * dr.x * swf->swf,
97 g * dr.y * swf->swf,
98 g * dr.z * swf->swf
99 };
100
101 efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x),
102 CVEC(pt_i->x), &force, NULL);
103 efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x),
104 CVEC(pt_j->x), &force, NULL);
105 efp_add_stress(&swf->dr, &force, &efp->stress);
106 }
107 return energy;
108 }
109
110 static double
disp_overlap(struct efp * efp,size_t fr_i_idx,size_t fr_j_idx,size_t pt_i_idx,size_t pt_j_idx,double s_ij,six_t ds_ij,double sum,const struct swf * swf)111 disp_overlap(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx,
112 size_t pt_i_idx, size_t pt_j_idx, double s_ij, six_t ds_ij,
113 double sum, const struct swf *swf)
114 {
115 const struct frag *fr_i = efp->frags + fr_i_idx;
116 const struct frag *fr_j = efp->frags + fr_j_idx;
117
118 const struct dynamic_polarizable_pt *pt_i =
119 fr_i->dynamic_polarizable_pts + pt_i_idx;
120 const struct dynamic_polarizable_pt *pt_j =
121 fr_j->dynamic_polarizable_pts + pt_j_idx;
122
123 vec_t dr = {
124 pt_j->x - pt_i->x - swf->cell.x,
125 pt_j->y - pt_i->y - swf->cell.y,
126 pt_j->z - pt_i->z - swf->cell.z
127 };
128
129 double r = vec_len(&dr);
130 double r2 = r * r;
131 double r6 = r2 * r2 * r2;
132
133 double ln_s = 0.0;
134 double damp = 1.0;
135
136 if (fabs(s_ij) > 1.0e-5) {
137 ln_s = log(fabs(s_ij));
138 damp = 1.0 - s_ij * s_ij * (1.0 - 2.0 * ln_s +
139 2.0 * ln_s * ln_s);
140 }
141
142 double energy = -4.0 / 3.0 * sum * damp / r6;
143
144 if (efp->do_gradient) {
145 vec_t force, torque_i, torque_j;
146
147 double t1 = -8.0 * sum / r6 / r2 * damp;
148 double t2 = -16.0 / 3.0 * sum / r6 * ln_s * ln_s * s_ij;
149
150 vec_t dr_i = vec_sub(CVEC(pt_i->x), CVEC(fr_i->x));
151 vec_t dr_j = vec_sub(CVEC(pt_j->x), CVEC(fr_j->x));
152
153 force.x = (t1 * dr.x - t2 * ds_ij.x) * swf->swf;
154 force.y = (t1 * dr.y - t2 * ds_ij.y) * swf->swf;
155 force.z = (t1 * dr.z - t2 * ds_ij.z) * swf->swf;
156
157 torque_i.x = swf->swf * (t1 * (dr.z * dr_i.y - dr.y * dr_i.z) +
158 t2 * ds_ij.a);
159 torque_i.y = swf->swf * (t1 * (dr.x * dr_i.z - dr.z * dr_i.x) +
160 t2 * ds_ij.b);
161 torque_i.z = swf->swf * (t1 * (dr.y * dr_i.x - dr.x * dr_i.y) +
162 t2 * ds_ij.c);
163
164 torque_j.x = swf->swf * (t1 * (dr.z * dr_j.y - dr.y * dr_j.z) +
165 t2 * (ds_ij.z * swf->dr.y - ds_ij.y * swf->dr.z) +
166 t2 * ds_ij.a);
167 torque_j.y = swf->swf * (t1 * (dr.x * dr_j.z - dr.z * dr_j.x) +
168 t2 * (ds_ij.x * swf->dr.z - ds_ij.z * swf->dr.x) +
169 t2 * ds_ij.b);
170 torque_j.z = swf->swf * (t1 * (dr.y * dr_j.x - dr.x * dr_j.y) +
171 t2 * (ds_ij.y * swf->dr.x - ds_ij.x * swf->dr.y) +
172 t2 * ds_ij.c);
173
174 six_atomic_add_xyz(efp->grad + fr_i_idx, &force);
175 six_atomic_add_abc(efp->grad + fr_i_idx, &torque_i);
176 six_atomic_sub_xyz(efp->grad + fr_j_idx, &force);
177 six_atomic_sub_abc(efp->grad + fr_j_idx, &torque_j);
178 efp_add_stress(&swf->dr, &force, &efp->stress);
179 }
180 return energy;
181 }
182
183 static double
disp_off(struct efp * efp,size_t fr_i_idx,size_t fr_j_idx,size_t pt_i_idx,size_t pt_j_idx,double sum,const struct swf * swf)184 disp_off(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, size_t pt_i_idx,
185 size_t pt_j_idx, double sum, const struct swf *swf)
186 {
187 const struct frag *fr_i = efp->frags + fr_i_idx;
188 const struct frag *fr_j = efp->frags + fr_j_idx;
189
190 const struct dynamic_polarizable_pt *pt_i =
191 fr_i->dynamic_polarizable_pts + pt_i_idx;
192 const struct dynamic_polarizable_pt *pt_j =
193 fr_j->dynamic_polarizable_pts + pt_j_idx;
194
195 vec_t dr = {
196 pt_j->x - pt_i->x - swf->cell.x,
197 pt_j->y - pt_i->y - swf->cell.y,
198 pt_j->z - pt_i->z - swf->cell.z
199 };
200
201 double r = vec_len(&dr);
202 double r2 = r * r;
203 double r6 = r2 * r2 * r2;
204
205 double energy = -4.0 / 3.0 * sum / r6;
206
207 if (efp->do_gradient) {
208 double r8 = r6 * r2;
209 double g = -8.0 * sum / r8;
210
211 vec_t force = {
212 g * dr.x * swf->swf,
213 g * dr.y * swf->swf,
214 g * dr.z * swf->swf
215 };
216
217 efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x),
218 CVEC(pt_i->x), &force, NULL);
219 efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x),
220 CVEC(pt_j->x), &force, NULL);
221 efp_add_stress(&swf->dr, &force, &efp->stress);
222 }
223 return energy;
224 }
225
226 static double
point_point_disp(struct efp * efp,size_t fr_i_idx,size_t fr_j_idx,size_t pt_i_idx,size_t pt_j_idx,double s,six_t ds,const struct swf * swf)227 point_point_disp(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx,
228 size_t pt_i_idx, size_t pt_j_idx, double s, six_t ds, const struct swf *swf)
229 {
230 struct frag *fr_i = efp->frags + fr_i_idx;
231 struct frag *fr_j = efp->frags + fr_j_idx;
232
233 const struct dynamic_polarizable_pt *pt_i =
234 fr_i->dynamic_polarizable_pts + pt_i_idx;
235 const struct dynamic_polarizable_pt *pt_j =
236 fr_j->dynamic_polarizable_pts + pt_j_idx;
237
238 double sum = 0.0;
239
240 for (size_t k = 0; k < ARRAY_SIZE(weights); k++) {
241 double tr_i = (pt_i->tensor[k].xx +
242 pt_i->tensor[k].yy +
243 pt_i->tensor[k].zz) / 3;
244 double tr_j = (pt_j->tensor[k].xx +
245 pt_j->tensor[k].yy +
246 pt_j->tensor[k].zz) / 3;
247 sum += weights[k] * tr_i * tr_j;
248 }
249
250 switch (efp->opts.disp_damp) {
251 case EFP_DISP_DAMP_TT:
252 return disp_tt(efp, fr_i_idx, fr_j_idx,
253 pt_i_idx, pt_j_idx, sum, swf);
254 case EFP_DISP_DAMP_OVERLAP:
255 return disp_overlap(efp, fr_i_idx, fr_j_idx,
256 pt_i_idx, pt_j_idx, s, ds, sum, swf);
257 case EFP_DISP_DAMP_OFF:
258 return disp_off(efp, fr_i_idx, fr_j_idx,
259 pt_i_idx, pt_j_idx, sum, swf);
260 }
261 assert(0);
262 }
263
264 /*
265 * Reference:
266 *
267 * Ivana Adamovic, Mark Gordon
268 *
269 * Dynamic polarizability, dispersion coefficient C6 and dispersion energy in
270 * the effective fragment potential method
271 *
272 * Mol. Phys. 103, 379 (2005)
273 */
274 double
efp_frag_frag_disp(struct efp * efp,size_t frag_i,size_t frag_j,const double * s,const six_t * ds)275 efp_frag_frag_disp(struct efp *efp, size_t frag_i, size_t frag_j,
276 const double *s, const six_t *ds)
277 {
278 double energy = 0.0;
279
280 struct frag *fr_i = efp->frags + frag_i;
281 struct frag *fr_j = efp->frags + frag_j;
282
283 size_t n_disp_i = fr_i->n_dynamic_polarizable_pts;
284 size_t n_disp_j = fr_j->n_dynamic_polarizable_pts;
285
286 struct swf swf = efp_make_swf(efp, fr_i, fr_j);
287
288 for (size_t ii = 0, idx = 0; ii < n_disp_i; ii++)
289 for (size_t jj = 0; jj < n_disp_j; jj++, idx++)
290 energy += point_point_disp(efp, frag_i, frag_j, ii, jj,
291 s[idx], ds[idx], &swf);
292
293 vec_t force = {
294 swf.dswf.x * energy,
295 swf.dswf.y * energy,
296 swf.dswf.z * energy
297 };
298
299 six_atomic_add_xyz(efp->grad + frag_i, &force);
300 six_atomic_sub_xyz(efp->grad + frag_j, &force);
301 efp_add_stress(&swf.dr, &force, &efp->stress);
302
303 return energy * swf.swf;
304 }
305
306 void
efp_update_disp(struct frag * frag)307 efp_update_disp(struct frag *frag)
308 {
309 for (size_t i = 0; i < frag->n_dynamic_polarizable_pts; i++) {
310 const struct dynamic_polarizable_pt *pt_in =
311 frag->lib->dynamic_polarizable_pts + i;
312 struct dynamic_polarizable_pt *pt_out =
313 frag->dynamic_polarizable_pts + i;
314
315 efp_move_pt(CVEC(frag->x), &frag->rotmat,
316 CVEC(pt_in->x), VEC(pt_out->x));
317 for (size_t j = 0; j < 12; j++) {
318 const mat_t *in = pt_in->tensor + j;
319 mat_t *out = pt_out->tensor + j;
320
321 efp_rotate_t2(&frag->rotmat, (const double *)in,
322 (double *)out);
323 }
324 }
325 }
326