lennard-jones-cuda/src/tasks/semiGCEfluctuations/semiGCEfluctuations.cpp at master · vlvovch/lennard-jones-cuda · GitHub

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include "MDSystem.h"
#include "NumberStatistics.h"

using namespace std;
using namespace SampleMoments;

int GetNsub(const MDSystem& syst, double fraction = 0.5) {
  double L = syst.L;
  int N = syst.m_config.N;

  double tL = L * pow(fraction, 1. / 3.);

  int ret = 0;

  for (int i = 0; i < 4 * N; i += 4) {
    double dx = (syst.h_Pos[i] - 0.5 * L);
    double dy = (syst.h_Pos[i + 1] - 0.5 * L);
    double dz = (syst.h_Pos[i + 2] - 0.5 * L);

    if (abs(dx) < 0.5 * tL && abs(dy) < 0.5 * tL && abs(dz) < 0.5 * tL)
      ret++;
  }

  return ret;
}

int main(int argc, char *argv[])
{
  int nev = 10000;
  vector<double> fracs;
  for (double tfr = 0.05; tfr <= 1.; tfr += 0.05)
    fracs.push_back(tfr);

  vector<NumberStatistics> StatsFlucs(fracs.size());

  int N = 512;// 4096;
  double Tst = 2.0;
  double rhost = 0.050;
  Tst = 1.312;
  rhost = 0.316;

  double dt = 0.005;
  int iterspreeq = 10000;// / 2;
  int itersstep = 200;// / 2;

  MDSystem::MDSystemConfiguration config;
  config.N = N;
  config.T0 = Tst;
  config.rho = rhost;
  config.useCUDA = true;

  MDSystem syst(config);

  syst.Reinitialize(config);


  // Equilibration phase
  syst.m_config.canonical = true;
  for (int i = 0; i < iterspreeq; ++i) {
    syst.Integrate(dt);
    //syst.RenormalizeVelocities();
  }


  // Production phase
  syst.m_config.canonical = false;
  for (int iN = 1; iN <= nev; iN++) {
    //syst.Reinitialize(config);

    //for (int i = 0; i < iterspreeq; ++i) {
    //  syst.Integrate(dt);
    //  syst.RenormalizeVelocities();
    //}

	//syst.RenormalizeVelocities();

    //syst.m_config.canonical = true;

    for (int i = 0; i < itersstep; ++i) {
      syst.Integrate(dt);
    }

    for (int i = 0; i < fracs.size(); ++i) {
      StatsFlucs[i].AddObservation(static_cast<double>(GetNsub(syst, fracs[i])));
    }

    //printf("%d ", iN);
    if (iN % 10 == 0) {
      for (int i = 0; i < fracs.size(); ++i) {
        printf("%15d %10lf +- %-10lf %10lf +- %-10lf %10lf +- %-10lf %10lf +- %-10lf\n",
          iN,
          StatsFlucs[i].GetMean(),
          StatsFlucs[i].GetMeanError(),
          StatsFlucs[i].GetScaledVariance(),
          StatsFlucs[i].GetScaledVarianceError(),
          StatsFlucs[i].GetScaledVariance() / (1. - fracs[i]),
          StatsFlucs[i].GetScaledVarianceError() / (1. - fracs[i]),
		      StatsFlucs[i].GetSkewness() / (1. - 2. * fracs[i]),
		      StatsFlucs[i].GetSkewnessError() / (1. - 2. * fracs[i])
          );
      }
      printf("\n");
    }
  }
  return 0;
}