blob: b1237c4c741929e211aaac1b536ef430dd283383 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
#include <Eigen/Core>
namespace Eigen {
template<typename Scalar> struct pow12_op EIGEN_EMPTY_STRUCT {
inline const Scalar operator() (const Scalar& a) const
{
Scalar b = a*a*a;
Scalar c = b*b;
return c*c;
}
template<typename PacketScalar>
inline const PacketScalar packetOp(const PacketScalar& a) const
{
PacketScalar b = ei_pmul(a, ei_pmul(a, a));
PacketScalar c = ei_pmul(b, b);
return ei_pmul(c, c);
}
};
template<typename Scalar>
struct ei_functor_traits<pow12_op<Scalar> >
{
enum {
Cost = 4 * NumTraits<Scalar>::MulCost,
PacketAccess = int(ei_packet_traits<Scalar>::size) > 1
};
};
} // namespace Eigen
using Eigen::pow12_op;
USING_PART_OF_NAMESPACE_EIGEN
#ifndef SCALAR
#define SCALAR float
#endif
#ifndef SIZE
#define SIZE 10000
#endif
#ifndef REPEAT
#define REPEAT 10000
#endif
typedef Matrix<SCALAR, Eigen::Dynamic, 1> Vec;
using namespace std;
SCALAR E_VDW(const Vec &interactions1, const Vec &interactions2)
{
return interactions2
.cwiseQuotient(interactions1)
.cwise(pow12_op<SCALAR>())
.sum();
}
int main()
{
//
// 1 2 3 4 ... (interactions)
// ka . . . . ...
// rab . . . . ...
// energy . . . . ...
// ... ... ... ... ... ...
// (variables
// for
// interaction)
//
Vec interactions1(SIZE), interactions2(SIZE); // SIZE is the number of vdw interactions in our system
// SetupCalculations()
SCALAR rab = 1.0;
interactions1.setConstant(2.4);
interactions2.setConstant(rab);
// Energy()
SCALAR energy = 0.0;
for (unsigned int i = 0; i<REPEAT; ++i) {
energy += E_VDW(interactions1, interactions2);
energy *= 1 + 1e-20 * i; // prevent compiler from optimizing the loop
}
cout << "energy = " << energy << endl;
}
|
