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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

/********************************************\
  *
  *  Sire - Molecular Simulation Framework
  *
  *  Copyright (C) 2006  Christopher Woods
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation; either version 2 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  *
  *  For full details of the license please see the COPYING file
  *  that should have come with this distribution.
  *
  *  You can contact the authors via the developer's mailing list
  *  at http://siremol.org
  *
\*********************************************/

#ifndef SIREVOL_SPACE_H
#define SIREVOL_SPACE_H

#include <QList>
#include <QVector>
#include <boost/tuple/tuple.hpp>

#include "SireBase/pairmatrix.hpp"

#include "SireUnits/dimensions.h"

#include "SireMaths/vector.h"
#include "SireMaths/distvector.h"

#include "Siren/object.h"
#include "Siren/objptr.hpp"

#include "coordgroup.h"

SIRE_BEGIN_HEADER

namespace SireMaths
{
class RanGenerator;
}

namespace SireVol
{

using SireMaths::Vector;
using SireMaths::DistVector;
using SireMaths::RanGenerator;

/** Define a distance matrix as being a PairMatrix of doubles */
typedef SireBase::PairMatrix<double> DistMatrix;

/** Define a vector field matrix as being a PairMatrix of DistVectors
    (these are Vectors that have the magnitude and direction separated out) */
typedef SireBase::PairMatrix<DistVector> DistVectorMatrix;

class Space;
typedef Siren::ObjPtr<Space> SpacePtr;

/** This pure virtual base class represents the volume of space within which a SimSystem
    resides. This class provides functions that calculate the distances between CoordGroups
    in this volume, and has most of the optimisation (since most of the hard-work
    double distance loops occur in this class!). Key overloaded classes that inherit
    Space are Cartesian, which represents infinite 3D cartesian space, and
    PeriodicBox which represents a 3D periodic box.

    As this class is used right in the heart of the double loop it is very important
    that it is not implemented in a way that requires a lot of virtual function calls.
    This means that the class is implemented to calculate the distances of all pairs of points
    between two CoordGroups in a single virtual function call. The results are held
    in a special matrix that can be queried via read-only, inline, non-virtual
    functions, with the items ordered in memory in the same order that you would get
    by iterating over all pairs (group1 in outer loop, group2 in inner loop). An advantage
    of this approach is that as the distances are calculated in one go, it is possible for
    the Space class to find out whether two CoordGroups are within the non-bonded cutoff
    distance before any further calculation is performed.

    To prevent continual reallocation, the Space class works on an already-allocated
    matrix class. This is only reallocated if it is found that it is not sufficiently
    large to hold the pair of CoordGroups.

    As a further optimisation, the distance matrix may be populated with the interpoint
    distances, or the square of the interatomic distances, or the 1 / distances
    or 1 / distances^2

    The inheritors of this class should be the only parts of this code where
    distance calculations are calculated between and within CoordGroups. This will
    allow you to change the space of the system here, and have that space used in the
    rest of the code that uses CoordGroups.

    This is a virtual class that is designed to be used with SharedPolyPointer.

    @author Christopher Woods
*/
class SIREVOL_EXPORT Space : public Siren::Extends<Space,Siren::Object>
{
public:
    Space();
    Space(const Space &other);

    ~Space();

    ///////////////////////////
    // Extends Siren::Object //
    ///////////////////////////
    
    static QString typeName();
    
    void stream(Siren::Stream &s);
    
    ////////////////////
    // SireVol::Space //
    ////////////////////

    /** Return the volume of the central box of this space. This
        throws an exception if it is not possible to calculate the
        volume of this space (e.g. it is an infinite space!) */
    virtual SireUnits::Dimension::Volume volume() const=0;

    /** Return a copy of this space with the volume of set to 'volume'
        - this will scale the space uniformly, keeping the center at
        the same location, to achieve this volume */
    virtual SpacePtr setVolume(SireUnits::Dimension::Volume volume) const=0;

    /** Change the volume of this space by 'delta' */
    SpacePtr changeVolume(SireUnits::Dimension::Volume delta) const;

    /** Calculate the distance between two points */
    virtual double calcDist(const Vector &point0, const Vector &point1) const=0;

    /** Calculate the distance squared between two points */
    virtual double calcDist2(const Vector &point0, const Vector &point1) const=0;

    /** Populate the matrix 'mat' with the distances between all of the
        points within a CoordGroup. This creates a symmetrical matrix,
        with a 0 diagonal. This returns the shortest distance between
        two points within the group. */
    virtual double calcDist(const CoordGroup &group, DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the distances^2 between all of the
        points within a CoordGroup. This creates a symmetrical matrix,
        with a 0 diagonal. This returns the shortest distance^2 between
        two points within the group. */
    virtual double calcDist2(const CoordGroup &group, DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the inverse distances between all of the
        points within a CoordGroup. This creates a symmetrical matrix,
        with a 0 diagonal. This returns the the largest inverse distance between
        two points within the group. */
    virtual double calcInvDist(const CoordGroup &group, DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the inverse distances^2 between all of the
        points within a CoordGroup. This creates a symmetrical matrix,
        with a 0 diagonal. This returns the the largest inverse distance^2 between
        two points within the group. */
    virtual double calcInvDist2(const CoordGroup &group, DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the distances between all of the
        points of the two CoordGroups. Return the shortest distance^2 between the two
        CoordGroups. */
    virtual double calcDist(const CoordGroup &group1, const CoordGroup &group2,
                            DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the distances^2 between all of the
        points of the two CoordGroups. Return the shortest distance^2 between the
        two CoordGroups. */
    virtual double calcDist2(const CoordGroup &group1, const CoordGroup &group2,
                             DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the inverse distances between all of the
        points of the two CoordGroups. Return the largest inverse distance between the two
        CoordGroups. */
    virtual double calcInvDist(const CoordGroup &group1, const CoordGroup &group2,
                               DistMatrix &distmat) const=0;

    /** Populate the matrix 'mat' with the inverse distances^2 between all of the
        points of the two CoordGroups. Return the largest inverse distance^2 between the two
        CoordGroups. */
    virtual double calcInvDist2(const CoordGroup &group1, const CoordGroup &group2,
                                DistMatrix &distmat) const=0;

    /** Calculate the distance vector between two points */
    virtual DistVector calcDistVector(const Vector &point0, 
                                      const Vector &point1) const=0;

    /** Populate the matrix 'distmat' with all of the interpoint distance vectors
        between all points within the CoordGroup. This is *not* a symmetrical matrix,
        as the direction from point A to point B is the negative of the 
        direction from point B to point A. This returns the shortest distance
        between two points in the group (that is not the self-self distance) */
    virtual double calcDistVectors(const CoordGroup &group,
                                   DistVectorMatrix &distmat) const=0;
                                       
    /** Populate the matrix 'distmat' between all the points of the two CoordGroups
        'group1' and 'group2' - the returned matrix has the vectors pointing
        from each point in 'group1' to each point in 'group2'. This returns
        the shortest distance between two points in the group */
    virtual double calcDistVectors(const CoordGroup &group1,
                                   const CoordGroup &group2,
                                   DistVectorMatrix &distmat) const=0;

    /** Calculate the angle between the passed three points. This should return
        the acute angle between the points, which should lie between 0 and 180 degrees */
    virtual SireUnits::Dimension::Angle calcAngle(const Vector &point0,
                                                  const Vector &point1,
                                                  const Vector &point2) const=0;

    /** Calculate the torsion angle between the passed four points. This should
        return the torsion angle measured clockwise when looking down the 
        torsion from point0-point1-point2-point3. This will lie between 0 and 360 
        degrees */
    virtual SireUnits::Dimension::Angle calcDihedral(const Vector &point0,
                                                     const Vector &point1,
                                                     const Vector &point2,
                                                     const Vector &point3) const=0;

    /** Return whether or not these two groups are definitely beyond the distance 'dist'.

        \warning Note 'beyond' does not mean definitely within the distance!
    */
    virtual bool beyond(double dist, const CoordGroup &group0,
                        const CoordGroup &group1) const=0;

    /** Return whether or not two groups that are enclosed by the AABoxes
        'aabox0' and 'aabox1' are beyond the cutoff distance 'dist'.
        
        \warning Note 'beyond' does not mean definitely within the distance!
    */
    virtual bool beyond(double dist, const AABox &aabox0, const AABox &aabox1) const=0;

    /** Return the minimum distance between the points in 'group0' and 'group1'.
        If this is a periodic space then this uses the minimum image convention
        (i.e. the minimum distance between the closest periodic replicas are
        used) */
    virtual double minimumDistance(const CoordGroup &group0,
                                   const CoordGroup &group1) const=0;

    /** Return the minimum distance between points within the group 'group'. */
    virtual double minimumDistance(const CoordGroup &group) const=0;

    /** Return whether or not this space is periodic */
    virtual bool isPeriodic() const=0;
    
    /** Return whether or not this space is cartesian
        ( x==y==z == 1 and all angle between x-y, y-z, x-z all 90 degrees ) */
    virtual bool isCartesian() const=0;

    /** Return the center of the box that contains the point 'p' assuming
        that the center for the central box is located at the origin */
    virtual Vector getBoxCenter(const Vector &p) const=0;

	/** Return the center of the box that contains the point 'p' assuming
        that the center for the central box is located at 'center' */
    virtual Vector getBoxCenter(const Vector &p, const Vector &center) const=0;

    /** Return a random point within this space, using the passed 
        random number generator to generate the necessary random numbers,
        and placing the center of the box at 'center' */
    virtual Vector getRandomPoint(const Vector &center, 
                                  const RanGenerator &generator) const=0;

    virtual Vector getRandomPoint(const Vector &center) const;

    virtual Vector getRandomPoint(const RanGenerator &generator) const;

    virtual Vector getRandomPoint() const;

    /** Return the minimum image copy of 'group' with respect to 'center'.
        For periodic spaces, this translates 'group' into the box that
        has its center at 'center' (i.e. returns the closest copy of
        'group' to 'center' according to the minimum image convention) */
    virtual CoordGroup getMinimumImage(const CoordGroup &group,
                                       const Vector &center) const=0;

    /** Return the minimum image copy of 'groups' with respect to 'center'.
        For periodic spaces, this translates 'groups' into the box that
        has its center at 'center' (i.e. returns the closest copy of
        each 'group' to 'center' according to the minimum image convention) */
    virtual CoordGroupArray getMinimumImage(const CoordGroupArray &groups,
                                            const Vector &center,
                                            bool translate_as_one=false) const=0;

    /** Return the minimum image copy of 'aabox' with respect to 'center'.  
        For periodic spaces, this returns the AABox translated into the 
        box that has its center at 'center' */
    virtual AABox getMinimumImage(const AABox &aabox, const Vector &center) const=0;
    
    /** Return the minimum image copy of 'point' with respect to 'center'. 
        For periodic spaces, this returns the point translated into the
        box that has its center at 'center' */
    virtual Vector getMinimumImage(const Vector &point, const Vector &center) const=0;

    /** Return a list of copies of CoordGroup 'group' that are within
        'distance' of the CoordGroup 'center', translating 'group' so that
        it has the right coordinates to be around 'center'. Note that multiple
        copies of 'group' may be returned if this is a periodic space and
        there are multiple periodic replicas of 'group' within 'dist' of
        'center'. The copies of 'group' are returned together with the
        minimum distance between that periodic replica and 'center'.

        If there are no periodic replicas of 'group' that are within
        'dist' of 'center', then an empty list is returned.
    */
    virtual QList< boost::tuple<double,CoordGroup> >
                    getCopiesWithin(const CoordGroup &group, const CoordGroup &center,
                                    double dist) const=0;

protected:
    Space& operator=(const Space &other);
    
    bool operator==(const Space &other) const;
    bool operator!=(const Space &other) const;
    
    void assertCompatible(const Space &other) const;
};

}

SIRE_EXPOSE_CLASS( SireVol::Space )

Q_DECLARE_METATYPE( SireVol::SpacePtr )

SIRE_EXPOSE_OBJECT_PTR( SireVol::SpacePtr, SireVol::Space )

SIRE_END_HEADER

#endif