source: src/samples/bspline.cpp@ 4571da

Last change on this file since 4571da was 716da7, checked in by Julian Iseringhausen <isering@…>, 14 years ago

Fix energy calculation.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1729 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 4.5 KB
Line 
1/**
2 * @file bspline.cpp
3 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
4 * @date Mon Nov 21 13:27:22 2011
5 *
6 * @brief B-Splines for molecular dynamics.
7 *
8 */
9
10#ifdef HAVE_CONFIG_H
11#include <config.h>
12#endif
13
14#ifndef BSPLINE_DEGREE
15#error BSPLINE_DEGREE not defined.
16#endif
17
18#if BSPLINE_DEGREE < 3 || BSPLINE_DEGREE > 6
19#error Please choose a B-Spline degree between 3 and 6
20#endif
21
22#include "base/helper.hpp"
23#include "base/math.hpp"
24#include "samples/bspline.hpp"
25
26#define POW(x,y) Helper::pow(x,y)
27
28using namespace VMG;
29
30BSpline::BSpline(const vmg_float& width) :
31 spline_nom((BSPLINE_DEGREE+1)/2),
32 spline_denom((BSPLINE_DEGREE+1)/2),
33 potential_nom((BSPLINE_DEGREE+1)/2+1),
34 potential_denom((BSPLINE_DEGREE+1)/2+1),
35 intervals((BSPLINE_DEGREE+1)/2),
36 R(width)
37{
38 for (unsigned int i=0; i<intervals.size(); ++i)
39 intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
40
41#if BSPLINE_DEGREE == 3
42
43 spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
44 spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
45
46 spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
47 spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));
48
49 potential_nom[0] = Polynomial(5,
50 0.0,
51 -195.0*POW(R,4),
52 0.0,
53 270.0*POW(R,2),
54 0.0,
55 -243.0);
56
57 potential_nom[1] = Polynomial(5,
58 2.0 * POW(R,5),
59 -405.0 * POW(R,4),
60 0.0,
61 810.0 * POW(R,2),
62 -810.0 * R,
63 243.0);
64
65 potential_nom[2] = Polynomial(0, -1.0);
66
67 potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
68 potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
69 potential_denom[2] = Polynomial(0, 1.0);
70
71 antid = 39.0 / (16.0 * R);
72
73#elif BSPLINE_DEGREE == 6
74
75 spline_nom[0] = Polynomial(5,
76 297.0 * POW(R,5),
77 0.0,
78 -2430.0 * POW(R,3),
79 0.0,
80 10935.0 * R,
81 -10935.0);
82
83 spline_nom[1] = Polynomial(5,
84 459.0 * POW(R,5),
85 2025.0 * POW(R,4),
86 -17010.0 * POW(R,3),
87 36450.0 * POW(R,2),
88 -32805.0 * R,
89 10935.0);
90
91 spline_nom[2] = Polynomial(5,
92 2187.0 * POW(R,5),
93 -10935.0 * POW(R,4),
94 21870.0 * POW(R,3),
95 -21870.0 * POW(R,2),
96 10935.0 * R,
97 -2187.0);
98
99 spline_denom[0] = Polynomial(0, 20.0 * Math::pi * POW(R,8));
100 spline_denom[1] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
101 spline_denom[2] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
102
103 potential_nom[0] = Polynomial(8,
104 0.0,
105 -956.0 * POW(R,7),
106 0.0,
107 2772.0 * POW(R,5),
108 0.0,
109 -6804.0 * POW(R,3),
110 0.0,
111 14580.0 * R,
112 -10935.0);
113
114 potential_nom[1] = Polynomial(8,
115 -5.0 * POW(R,8),
116 -5676.0 * POW(R,7),
117 0.0,
118 12852.0 * POW(R,5),
119 28350.0 * POW(R,4),
120 -142884.0 * POW(R,3),
121 204120.0 * POW(R,2),
122 -131220.0 * R,
123 32805.0);
124
125 potential_nom[2] = Polynomial(8,
126 169.0 * POW(R,8),
127 -2916.0 * POW(R,7),
128 0.0,
129 20412.0 * POW(R,5),
130 -51030.0 * POW(R,4),
131 61236.0 * POW(R,3),
132 -40824.0 * POW(R,2),
133 14580.0 * R,
134 -2187.0);
135
136 potential_nom[3] = Polynomial(0, -1.0);
137
138 potential_denom[0] = Polynomial(0, 280.0 * POW(R,8));
139 potential_denom[1] = Polynomial(0, 1680.0 * POW(R,8));
140 potential_denom[2] = Polynomial(0, 560.0 * POW(R,8));
141 potential_denom[3] = Polynomial(0, 0.0);
142
143 antid = 239.0 / (70.0 * R);
144
145#else
146#error B-Spline degree not supported. Choose 3 or 6.
147#endif /* BSPLINE_DEGREE */
148}
Note: See TracBrowser for help on using the repository browser.