1 /*
2 
3 Copyright (C) 2007 Lucas K. Wagner
4 
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9 
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 GNU General Public License for more details.
14 
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 
19 */
20 
21 
22 #ifndef MOLECULAR_SAMPLE_H_INCLUDED
23 #define MOLECULAR_SAMPLE_H_INCLUDED
24 
25 
26 #include "Sample_point.h"
27 class Wavefunction;
28 #include "Molecular_system.h"
29 class Sample_storage;
30 
31 
32 class Molecular_sample : public Sample_point
33 {
34 public:
35 
~Molecular_sample()36   ~Molecular_sample()
37   {}
38 
39 
40   void init(System * sys);
41 
42   void randomGuess();
43 
electronSize()44   int electronSize()
45   {
46     return nelectrons;
47   }
ionSize()48   int ionSize()
49   {
50     return parent->ions.size();
51   }
centerSize()52   int centerSize()
53   {
54     return ionSize();
55   }
getIonCharge(const int ion)56   doublevar getIonCharge(const int ion)
57   {
58     return parent->ions.charge(ion);
59   }
60 
61 
62   /*!
63     Moves electron e
64     position should be an array of length at least 3
65 
66   */
67   void setElectronPos(const int e,const Array1 <doublevar> & position);
68   virtual void setElectronPosNoNotify(const int e,
69                                     const Array1 <doublevar> & pos);
70 
71   void getElectronPos(const int e, Array1 <doublevar> & R);
72 
73   void moveIon(const int ion, const Array1 <doublevar> & r);
74 
getIonPos(const int ion,Array1<doublevar> & r)75   void getIonPos(const int ion, Array1 <doublevar> & r)
76   {
77     assert(r.GetDim(0) >= 3);
78 
79     for(int i=0; i<3; i++)
80     {
81       r(i)=parent->ions.r(i,ion);
82     }
83   }
84 
85 
86 
87   void updateEIDist();
88   void updateEEDist();
updateECDist()89   void updateECDist()
90   {
91     updateEIDist();
92   }
93 
94 
95 
96 
getECDist(const int e,const int cent,Array1<doublevar> & distance)97   void getECDist(const int e, const int cent,
98                  Array1 <doublevar> & distance)
99   {
100     getEIDist(e,cent, distance);
101   }
102 
getEIDist(const int e,const int ion,Array1<doublevar> & distance)103   void getEIDist(const int e,const int ion, Array1 <doublevar> & distance)
104   {
105     assert( distance.GetDim(0) >= 5 );
106     assert( ! ionDistStale(e));
107     for(int i=0; i<5; i++)
108     {
109       distance(i)=iondist(i,ion,e);
110     }
111   }
112 
getEIDist_temp(const int e,int ion,Array1<doublevar> & distance)113   virtual int getEIDist_temp(const int e, int ion,
114 			      Array1 <doublevar> & distance) {
115     assert(distance.GetDim(0)>=5);
116     distance(1)=0;
117     for(int d=0; d< 3; d++) {
118       distance(d+2)=elecpos(e,d)-parent->ions.r(d,ion);
119       distance(1)+=distance(d+2)*distance(d+2);
120     }
121     distance(0)=sqrt(distance(1));
122     return 1;
123   }
124 
getEEDist(const int e1,const int e2,Array1<doublevar> & distance)125   void getEEDist(const int e1,const int e2, Array1 <doublevar> & distance)
126   {
127     //cout << "getEEDist" << endl;
128     assert(distance.GetDim(0) >= 5);
129     assert( ! elecDistStale(e1));
130     assert( e1 < e2 );
131     for(int i=0; i< 5; i++)
132     {
133       distance(i)=pointdist(i,e1,e2);
134     }
135   }
136 
137   void rawOutput(ostream &);
138   void rawInput(istream &);
139 
140   void generateStorage(Sample_storage *& );
141   void saveUpdate(int, Sample_storage * );
142   void restoreUpdate(int, Sample_storage *);
143 
144 private:
145 
146   int nelectrons;
147 
148   Array2 <doublevar> elecpos; //electron positions
149 
150   Array3 <doublevar> iondist;
151   Array1 <int> elecDistStale;
152   Array1 <int> ionDistStale;
153   Array3 <doublevar> pointdist;  //this is an upper triangular matrix;
154                                  //the lower is currently wasted
155 
156   Molecular_system * parent;
157 
158 };
159 
160 
161 
162 #endif //MOLECULAR_SAMPLE_H_INCLUDED
163 //-------------------------------------------------------------------------
164