1#include <iostream>
2#include <Eigen/Core>
3
4using namespace Eigen;
5
6#ifndef SCALAR
7#define SCALAR float
8#endif
9
10#ifndef SIZE
11#define SIZE 10000
12#endif
13
14#ifndef REPEAT
15#define REPEAT 10000
16#endif
17
18typedef Matrix<SCALAR, Eigen::Dynamic, 1> Vec;
19
20using namespace std;
21
22SCALAR E_VDW(const Vec &interactions1, const Vec &interactions2)
23{
24  return (interactions2.cwise()/interactions1)
25         .cwise().cube()
26         .cwise().square()
27         .cwise().square()
28         .sum();
29}
30
31int main()
32{
33  //
34  //          1   2   3   4  ... (interactions)
35  // ka       .   .   .   .  ...
36  // rab      .   .   .   .  ...
37  // energy   .   .   .   .  ...
38  // ...     ... ... ... ... ...
39  // (variables
40  //    for
41  // interaction)
42  //
43  Vec interactions1(SIZE), interactions2(SIZE); // SIZE is the number of vdw interactions in our system
44  // SetupCalculations()
45  SCALAR rab = 1.0;
46  interactions1.setConstant(2.4);
47  interactions2.setConstant(rab);
48
49  // Energy()
50  SCALAR energy = 0.0;
51  for (unsigned int i = 0; i<REPEAT; ++i) {
52    energy += E_VDW(interactions1, interactions2);
53    energy *= 1 + 1e-20 * i; // prevent compiler from optimizing the loop
54  }
55  cout << "energy = " << energy << endl;
56}
57