6
6
7
7
namespace libcppe {
8
8
9
- Eigen::VectorXd BMatrix::apply_direct (const Eigen::VectorXd& induced_moments) {
10
- Eigen::VectorXd ret = Eigen::VectorXd::Zero (induced_moments.size ());
11
- #pragma omp parallel for
12
- for (int i = 0 ; i < m_n_polsites; ++i) {
13
- int l = i * 3 ;
14
- Potential& pot1 = m_polsites[i];
15
- double * retx = &ret[l + 0 ];
16
- double * rety = &ret[l + 1 ];
17
- double * retz = &ret[l + 2 ];
18
- for (int j = 0 ; j < m_n_polsites; ++j) {
19
- int m = j * 3 ;
20
- Potential& pot2 = m_polsites[j];
21
- if (pot1.excludes_site (pot2.index ) || i == j) continue ;
22
- Eigen::Vector3d diff = pot2.get_site_position () - pot1.get_site_position ();
23
- Eigen::VectorXd T2;
24
- if (m_options.damp_induced ) {
25
- Polarizability& alpha_i = pot1.get_polarizability ();
26
- Polarizability& alpha_j = pot2.get_polarizability ();
27
- double v = m_options.damping_factor_induced /
28
- std::pow (alpha_i.get_isotropic_value () * alpha_j.get_isotropic_value (),
29
- 1.0 / 6.0 );
30
- T2 = tensors::T2_damp_thole (diff, v);
31
- } else {
32
- T2 = tensors::T2 (diff);
33
- }
34
- double fx = induced_moments[m + 0 ];
35
- double fy = induced_moments[m + 1 ];
36
- double fz = induced_moments[m + 2 ];
37
- // inline matrix-vector product, faster than unfolding and multiplying
38
- *retx -= T2[0 ] * fx + T2[1 ] * fy + T2[2 ] * fz;
39
- *rety -= T2[1 ] * fx + T2[3 ] * fy + T2[4 ] * fz;
40
- *retz -= T2[2 ] * fx + T2[4 ] * fy + T2[5 ] * fz;
41
- }
42
- ret.segment <3 >(l) += m_alpha_inverse[i] * induced_moments.segment <3 >(l);
43
- }
44
- return ret;
45
- }
46
-
47
- Eigen::VectorXd BMatrix::apply_fast_summation (const Eigen::VectorXd& induced_moments,
48
- std::string scheme) {
9
+ Eigen::VectorXd BMatrix::apply (const Eigen::VectorXd& induced_moments) {
49
10
Eigen::VectorXd ret = Eigen::VectorXd::Zero (induced_moments.size ());
50
11
51
12
int n_crit = m_options.tree_ncrit ;
@@ -59,33 +20,27 @@ Eigen::VectorXd BMatrix::apply_fast_summation(const Eigen::VectorXd& induced_mom
59
20
S[i * 3 + 1 ] = s (1 );
60
21
S[i * 3 + 2 ] = s (2 );
61
22
}
62
- std::shared_ptr<Tree<1 , 3 >> tree = build_shared_tree<1 , 3 >(
63
- m_positions.data (), S.data (), m_n_polsites, n_crit, order, theta, m_exclusions);
23
+
24
+ double damping = 0.0 ;
25
+ if (m_options.damp_induced ) {
26
+ damping = m_options.damping_factor_induced ;
27
+ }
28
+ std::shared_ptr<Tree<1 , 3 >> tree =
29
+ build_shared_tree<1 , 3 >(m_polsites, S.data (), n_crit, order, theta, damping);
64
30
std::vector<double > induced_fields_v (3 * m_n_polsites);
31
+ auto scheme = m_options.summation_induced_fields ;
65
32
if (scheme == " fmm" ) {
66
33
tree->compute_field_fmm (induced_fields_v.data ());
67
34
} else {
68
- throw std::runtime_error ( " No such summation scheme " + scheme );
35
+ tree-> compute_field_exact (induced_fields_v. data () );
69
36
}
70
37
Eigen::VectorXd induced_fields =
71
38
Eigen::Map<Eigen::VectorXd>(induced_fields_v.data (), induced_fields_v.size ());
72
39
return induced_fields + apply_diagonal (induced_moments);
73
40
}
74
41
75
- Eigen::VectorXd BMatrix::apply (const Eigen::VectorXd& induced_moments) {
76
- auto scheme = m_options.summation_induced_fields ;
77
- if (scheme == " direct" ) {
78
- return apply_direct (induced_moments);
79
- } else if (scheme == " fmm" ) {
80
- return apply_fast_summation (induced_moments, scheme);
81
- } else {
82
- throw std::invalid_argument (" Invalid summation scheme for induced fields provided." );
83
- }
84
- }
85
-
86
42
Eigen::VectorXd BMatrix::apply_diagonal_inverse (const Eigen::VectorXd& in) {
87
43
Eigen::VectorXd ret = Eigen::VectorXd::Zero (3 * m_n_polsites);
88
-
89
44
#pragma omp parallel for
90
45
for (int i = 0 ; i < m_n_polsites; ++i) {
91
46
int l = i * 3 ;
@@ -106,45 +61,6 @@ Eigen::VectorXd BMatrix::apply_diagonal(const Eigen::VectorXd& in) {
106
61
return ret;
107
62
}
108
63
109
- Eigen::VectorXd BMatrix::gauss_seidel_update (Eigen::VectorXd induced_moments,
110
- const Eigen::VectorXd& total_fields) {
111
- #pragma omp parallel for
112
- for (int i = 0 ; i < m_n_polsites; ++i) {
113
- Eigen::Vector3d Ftmp = Eigen::Vector3d::Zero ();
114
- int l = i * 3 ;
115
- Potential& pot1 = m_polsites[i];
116
- for (int j = 0 ; j < m_n_polsites; ++j) {
117
- int m = 3 * j;
118
- Potential& pot2 = m_polsites[j];
119
- if (pot1.excludes_site (pot2.index ) || i == j) {
120
- continue ;
121
- }
122
- Eigen::Vector3d diff = pot2.get_site_position () - pot1.get_site_position ();
123
- Eigen::VectorXd T2;
124
- if (m_options.damp_induced ) {
125
- Polarizability& alpha_i = pot1.get_polarizability ();
126
- Polarizability& alpha_j = pot2.get_polarizability ();
127
- double v = m_options.damping_factor_induced /
128
- std::pow (alpha_i.get_isotropic_value () * alpha_j.get_isotropic_value (),
129
- 1.0 / 6.0 );
130
- T2 = tensors::T2_damp_thole (diff, v);
131
- } else {
132
- T2 = tensors::T2 (diff);
133
- }
134
- double fx = induced_moments[m + 0 ];
135
- double fy = induced_moments[m + 1 ];
136
- double fz = induced_moments[m + 2 ];
137
- // inline matrix-vector product, faster than unfolding and multiplying
138
- Ftmp[0 ] += T2[0 ] * fx + T2[1 ] * fy + T2[2 ] * fz;
139
- Ftmp[1 ] += T2[1 ] * fx + T2[3 ] * fy + T2[4 ] * fz;
140
- Ftmp[2 ] += T2[2 ] * fx + T2[4 ] * fy + T2[5 ] * fz;
141
- }
142
- Ftmp += total_fields.segment <3 >(l);
143
- induced_moments.segment <3 >(l) = pot1.get_polarizability ().get_matrix () * Ftmp;
144
- }
145
- return induced_moments;
146
- }
147
-
148
64
Eigen::MatrixXd BMatrix::to_dense_matrix () {
149
65
Eigen::MatrixXd B = Eigen::MatrixXd::Zero (m_n_polsites * 3 , m_n_polsites * 3 );
150
66
#pragma omp parallel for
0 commit comments