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