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