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 "balance.h"
28 #include "private.h"
29
30 static const double quad_fact[12] = {
31 0.72086099022968040154e-02, 0.17697067815034886394e-01,
32 0.30660908596251749739e-01, 0.48381293256249884995e-01,
33 0.74878830420650517080e-01, 0.11806515901361630228e+00,
34 0.19535413832209084204e+00, 0.35055692324483221824e+00,
35 0.71577113554429568336e+00, 1.81409759976323969729e+00,
36 6.97923445114870823247e+00, 8.32480938829658453917e+01
37 };
38
39 static const double quad_freq[12] = {
40 0.77932702233253671082e-05, 0.22821071773724297874e-03,
41 0.15211319247778075068e-02, 0.60833919905855461032e-02,
42 0.19223967039304198946e-01, 0.54392829363594207533e-01,
43 0.14891668800412559598e+00, 0.42134903703482291779e+00,
44 1.33149401066630080343e+00, 5.32498192172462030801e+00,
45 3.54935126637048206534e+01, 1.03935828835455831714e+03
46 };
47
48 static double
quadrature(const mat_t * tensor,size_t i,size_t j,double de)49 quadrature(const mat_t *tensor, size_t i, size_t j, double de)
50 {
51 double m, sum = 0.0;
52
53 for (int k = 0; k < 12; k++) {
54 m = mat_get(tensor + k, i, j);
55 sum += m * quad_fact[k] / (de * de + quad_freq[k]);
56 }
57 return sum * de;
58 }
59
60 static double
get_dip_int(struct efp * efp,size_t i_occ,size_t i_vir,size_t axis)61 get_dip_int(struct efp *efp, size_t i_occ, size_t i_vir, size_t axis)
62 {
63 size_t idx, size;
64
65 size = efp->n_ai_core + efp->n_ai_act + efp->n_ai_vir;
66 idx = axis * size * size + i_occ * size +
67 (efp->n_ai_core + efp->n_ai_act) + i_vir;
68
69 return efp->ai_dipole_integrals[idx];
70 }
71
72 static double
compute_ai_disp_pt(struct efp * efp,size_t fr_idx,size_t pt_idx)73 compute_ai_disp_pt(struct efp *efp, size_t fr_idx, size_t pt_idx)
74 {
75 struct frag *frag;
76 struct dynamic_polarizable_pt *pt;
77 double sum = 0.0;
78 size_t ncoreact;
79
80 frag = efp->frags + fr_idx;
81 pt = frag->dynamic_polarizable_pts + pt_idx;
82 ncoreact = efp->n_ai_core + efp->n_ai_act;
83
84 for (size_t i_vir = 0; i_vir < efp->n_ai_vir; i_vir++) {
85 double e_vir = efp->ai_orbital_energies[ncoreact + i_vir];
86
87 for (size_t i_occ = 0; i_occ < ncoreact; i_occ++) {
88 double e_occ = efp->ai_orbital_energies[i_occ];
89
90 for (size_t i = 0; i < 3; i++) {
91 double di = get_dip_int(efp,
92 i_occ, i_vir, i);
93
94 for (size_t j = 0; j < 3; j++) {
95 double dj = get_dip_int(efp,
96 i_occ, i_vir, j);
97
98 sum += di * dj * quadrature(pt->tensor,
99 i, j, e_vir - e_occ);
100 }
101 }
102 }
103 }
104 return -sum / PI;
105 }
106
107 static void
compute_ai_disp_range(struct efp * efp,size_t from,size_t to,void * data)108 compute_ai_disp_range(struct efp *efp, size_t from, size_t to, void *data)
109 {
110 double energy = 0.0;
111
112 (void)data;
113
114 #ifdef _OPENMP
115 #pragma omp parallel for schedule(dynamic) reduction(+:energy)
116 #endif
117 for (size_t i = from; i < to; i++) {
118 size_t n_pt = efp->frags[i].n_dynamic_polarizable_pts;
119
120 for (size_t j = 0; j < n_pt; j++)
121 energy += compute_ai_disp_pt(efp, i, j);
122 }
123 efp->energy.ai_dispersion += energy;
124 }
125
126 enum efp_result
efp_compute_ai_disp(struct efp * efp)127 efp_compute_ai_disp(struct efp *efp)
128 {
129 if (!(efp->opts.terms & EFP_TERM_AI_DISP))
130 return EFP_RESULT_SUCCESS;
131
132 if (efp->do_gradient) {
133 efp_log("gradient for AI/EFP dispersion is not implemented");
134 return EFP_RESULT_FATAL;
135 }
136
137 efp_balance_work(efp, compute_ai_disp_range, NULL);
138 efp_allreduce(&efp->energy.ai_dispersion, 1);
139
140 return EFP_RESULT_SUCCESS;
141 }
142