1 #include <octave/oct.h>
2 #include <octave/parse.h>
3
min_image(Matrix P,double L)4 Matrix min_image (Matrix P, double L)
5 {
6 // Returns the image of a particle inside a square box of side L centered in the origin
7 Matrix Pn = P;
8 for (octave_idx_type ip = 0; ip < P.rows (); ip++)
9 for (octave_idx_type dim = 0; dim < P.columns (); dim++)
10 {
11 if (P(ip,dim) > L/2)
12 Pn.fill(P(ip,dim) - L * labs (lrint (P(ip,dim) / L)), ip, dim, ip, dim);
13
14 else if (P(ip,dim) < -L/2)
15 Pn.fill(P(ip,dim) + L * labs (lrint (P(ip,dim) / L)), ip, dim, ip, dim);
16 }
17 return Pn;
18 }
19
20 /* When this version is used, the interaction function must depend on the
21 relative position and velocities only.
22 */
23 DEFUN_DLD (verletstep_boxed, args, nargout, "Verlet velocity step in periodic space")
24 {
25 int nargin = args.length();
26 octave_value_list retval;
27
28 unsigned int fcn_str = 0;
29
30 if (nargin < 6)
31 print_usage ();
32 else
33 {
34 // Allocate inputs
35 octave_function *fcn;
36 if (args(4).is_function_handle () || args(4).is_inline_function ())
37 fcn_str = 0;
38 else if (args(4).is_string ())
39 fcn_str = 1;
40 else
41 error ("verletstep: expected string,"," inline or function handle");
42
43 Matrix P = args(0).matrix_value ();
44 Matrix V = args(1).matrix_value ();
45 Matrix M = args(2).matrix_value ();
46 double dt = args(3).scalar_value ();
47
48 dim_vector dv = P.dims();
49 octave_idx_type Nparticles = dv(0);
50 octave_idx_type Ndim = dv(1);
51
52 octave_value_list newargs;
53 Matrix A (dim_vector(Nparticles, Ndim),0);
54
55 // Box containing the particles
56 double L = args(5).scalar_value ();
57
58 // Evaluate interaction function
59 for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
60 {
61
62 for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
63 {
64 newargs(0) = min_image (P.row(ipart) - P.row(jpart), L);
65 newargs(1) = V.row(ipart) - V.row(jpart);
66
67 if (fcn_str)
68 retval = octave::feval (args (4).string_value (), newargs, nargout);
69 else
70 {
71 fcn = args(4).function_value ();
72 if (! error_state)
73 retval = octave::feval (fcn, newargs, nargout);
74 }
75
76 A.insert (retval(0).row_vector_value () +
77 A.row (ipart), ipart, 0);
78 A.insert (retval(1).row_vector_value () +
79 A.row (jpart), jpart, 0);
80
81 }
82 }
83
84 for (octave_idx_type i = 0; i < A.rows (); i++)
85 A.insert (A.row (i) / M(i), i, 0);
86
87 // Integrate by half time step velocity and full step position
88 V = V + A * dt/2;
89 P = P + V * dt;
90
91 A.fill(0);
92
93 // Evaluate interaction function
94 // Evaluate interaction function
95 for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
96 {
97
98 for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
99 {
100 newargs(0) = min_image (P.row(ipart) - P.row(jpart), L);
101 newargs(1) = V.row(ipart) - V.row(jpart);
102
103 if (fcn_str)
104 retval = octave::feval (args (4).string_value (), newargs, nargout);
105 else
106 {
107 fcn = args(4).function_value ();
108 if (! error_state)
109 retval = octave::feval (fcn, newargs, nargout);
110 }
111
112 A.insert (retval(0).row_vector_value () +
113 A.row (ipart), ipart, 0);
114 A.insert (retval(1).row_vector_value () +
115 A.row (jpart), jpart, 0);
116
117 }
118 }
119 for (octave_idx_type i = 0; i < A.rows (); i++)
120 A.insert (A.row (i) / M(i), i, 0);
121
122 // Integrate velocity for the rest of the time step
123 V = V + A * dt/2;
124
125 // Put all particles in a box of side L
126 P = min_image (P, L);
127
128 retval(0) = P;
129 retval(1) = V;
130 retval(2) = A;
131 }
132
133 return retval;
134 }
135
136