ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Atom.hpp
(Generate patch)

Comparing:
branches/mmeineke/OOPSE/libmdtools/Atom.hpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Atom.hpp (file contents), Revision 999 by chrisfen, Fri Jan 30 15:01:09 2004 UTC

# Line 1 | Line 1
1   #ifndef _ATOM_H_
2   #define _ATOM_H_
3  
4 < #include <cstring>
5 < #include <cstdlib>
4 > #include <string.h>
5 > #include <stdlib.h>
6   #include <iostream>
7  
8 + #include "SimState.hpp"
9 +
10   class Atom{
11   public:
10  Atom(int theIndex) {
11    c_n_hyd = 0;
12    has_dipole = 0;
13    is_VDW = 0;
14    is_LJ = 0;
12  
13 <    index = theIndex;
14 <    offset = 3 * index;
18 <    offsetX = offset;
19 <    offsetY = offset+1;
20 <    offsetZ = offset+2;
13 >  Atom(int theIndex, SimState* theConfig );
14 >  virtual ~Atom() {}
15  
16 <    Axx = index*9;
23 <    Axy = Axx+1;
24 <    Axz = Axx+2;
25 <    
26 <    Ayx = Axx+3;
27 <    Ayy = Axx+4;
28 <    Ayz = Axx+5;
16 >  virtual void setCoords(void);
17  
18 <    Azx = Axx+6;
19 <    Azy = Axx+7;
32 <    Azz = Axx+8;
33 <  }
34 <  virtual ~Atom() {}
18 >  void getPos( double theP[3] );
19 >  void setPos( double theP[3] );
20  
21 <  static void createArrays (int nElements) {
22 <    int i;
38 <    
39 <    pos = new double[nElements*3];
40 <    vel = new double[nElements*3];
41 <    frc = new double[nElements*3];
42 <    trq = new double[nElements*3];
43 <    Amat = new double[nElements*9];
44 <    mu = new double[nElements];
45 <    ul = new double[nElements*3];
21 >  void getVel( double theV[3] );
22 >  void setVel( double theV[3] );
23  
24 <    // init directional values to zero
25 <    
49 <    for( i=0; i<nElements; i++){
50 <      trq[i] = 0.0;
51 <      trq[i+1] = 0.0;
52 <      trq[i+2] = 0.0;
53 <      
54 <      Amat[i] = 1.0;
55 <      Amat[i+1] = 0.0;
56 <      Amat[i+2] = 0.0;
57 <      
58 <      Amat[i+3] = 0.0;
59 <      Amat[i+4] = 1.0;
60 <      Amat[i+5] = 0.0;
61 <      
62 <      Amat[i+6] = 0.0;
63 <      Amat[i+7] = 0.0;
64 <      Amat[i+8] = 1.0;
65 <      
66 <      mu[i] = 0.0;    
67 <      
68 <      ul[i] = 0.0;
69 <      ul[i+1] = 0.0;
70 <      ul[i+2] = 0.0;
71 <    }
72 <  }
73 <  static void destroyArrays(void) {
74 <    delete[] pos;
75 <    delete[] vel;
76 <    delete[] frc;
77 <    delete[] trq;
78 <    delete[] Amat;
79 <    delete[] mu;
80 <  }
24 >  void getFrc( double theF[3] );
25 >  void addFrc( double theF[3] );
26  
82  static double* getPosArray( void ) { return pos; }
83  static double* getVelArray( void ) { return vel; }
84  static double* getFrcArray( void ) { return frc; }
85  static double* getTrqArray( void ) { return trq; }
86  static double* getAmatArray( void ) { return Amat; }
87  static double* getMuArray( void ) { return mu; }
88  static double* getUlArray( void ) { return ul; }
89  
90  double getX() const {return pos[offsetX];}
91  double getY() const {return pos[offsetY];}
92  double getZ() const {return pos[offsetZ];}
93  void setX(double x) {pos[offsetX] = x;}
94  void setY(double y) {pos[offsetY] = y;}
95  void setZ(double z) {pos[offsetZ] = z;}
96  
97  double get_vx() const  {return vel[offsetX];}
98  double get_vy() const  {return vel[offsetY];}
99  double get_vz() const  {return vel[offsetZ];}
100  void set_vx(double vx) {vel[offsetX] = vx;}
101  void set_vy(double vy) {vel[offsetY] = vy;}
102  void set_vz(double vz) {vel[offsetZ] = vz;}
103  
104  double getFx() const   {return frc[offsetX];}
105  double getFy() const   {return frc[offsetY];}
106  double getFz() const   {return frc[offsetZ];}
107  void addFx(double add) {frc[offsetX] += add;}
108  void addFy(double add) {frc[offsetY] += add;}
109  void addFz(double add) {frc[offsetZ] += add;}
27    virtual void zeroForces() = 0;
28  
29    double getMass() const {return c_mass;}
30    void setMass(double mass) {c_mass = mass;}
31 +
32 +  double getEamRcut() const {return myEamRcut;}
33 +  void setEamRcut(double eamRcut) {myEamRcut = eamRcut;}
34    
35    double getSigma() const {return c_sigma;}
36    void setSigma(double sigma) {c_sigma = sigma;}
# Line 122 | Line 42 | class Atom{ (public)
42    void setCovalent(double covalent) {c_covalent = covalent;}
43    
44    int getIndex() const {return index;}
45 <  void setIndex(int theIndex) {
126 <    index = theIndex;
127 <    offset = index*3;
128 <    offsetX = offset;
129 <    offsetY = offset+1;
130 <    offsetZ = offset+2;
131 <
132 <    Axx = index*9;
133 <    Axy = Axx+1;
134 <    Axz = Axx+2;
135 <    
136 <    Ayx = Axx+3;
137 <    Ayy = Axx+4;
138 <    Ayz = Axx+5;
139 <
140 <    Azx = Axx+6;
141 <    Azy = Axx+7;
142 <    Azz = Axx+8;
143 <  }
144 <
45 >  void setIndex(int theIndex);
46    char *getType() {return c_name;}
47    void setType(char * name) {strcpy(c_name,name);}
48    
# Line 162 | Line 63 | class Atom{ (public)
63    void setLJ( void )        { is_LJ = 1; is_VDW = 0; }
64    int isLJ( void )    { return is_LJ; }
65  
66 <  void seVDW( void )        { is_VDW = 1; is_LJ = 0; }
66 >  void setVDW( void )        { is_VDW = 1; is_LJ = 0; }
67    int isVDW( void )    { return is_VDW; }
68 +
69 +  void setEAM( void ) { is_EAM = 1; }
70 +  int  isEAM( void ) { return is_EAM; }
71  
72 +  void setCharged( void ) { is_charged = 1; }
73 +  int isCharged( void ) { return is_charged; }
74 +
75    virtual int isDirectional( void ) = 0;
76  
170  static double* pos; // the position array
171  static double* vel; // the velocity array
172  static double* frc; // the forc array
173  static double* trq; // the torque vector  ( space fixed )
174  static double* Amat; // the rotation matrix
175  static double* mu; // the dipole moment array
176  static double* ul; // the lab frame unit directional vector
77  
78   protected:
79    
80 +  SimState* myConfig;
81 +
82 +  double* pos; // the position array
83 +  double* vel; // the velocity array
84 +  double* frc; // the forc array
85 +  double* trq; // the torque vector  ( space fixed )
86 +  double* Amat; // the rotation matrix
87 +  double* mu;   // the array of dipole moments
88 +  double* ul;   // the lab frame unit directional vector
89 +
90    double c_mass; /* the mass of the atom in amu */
91    double c_sigma; /* the sigma parameter for van der walls interactions */
92    double c_epslon; /* the esplon parameter for VDW interactions */
93    double c_covalent; // The covalent radius of the atom.
94  
95 +  double myEamRcut; // Atom rcut for eam defined by the forcefield.
96 +
97    int index; /* set the atom's index */
98    int offset; // the atom's offset in the storage array
99    int offsetX, offsetY, offsetZ;
# Line 196 | Line 108 | class Atom{ (public)
108    int c_n_hyd; // the number of hydrogens bonded to the atom
109    
110    int has_dipole; // dipole boolean
111 <  int is_VDW;    // VDW boolean
112 <  int is_LJ;    // LJ boolean
111 >  int is_VDW;     // VDW boolean
112 >  int is_LJ;      // LJ boolean
113 >  int is_EAM;     // EAM boolean
114 >  int is_charged;   // isCharge boolean
115  
116 +  bool hasCoords;
117 +
118   #ifdef IS_MPI
119    int myGlobalIndex;
120   #endif
121    
122   };
123  
208
209
124   class GeneralAtom : public Atom{
125  
126   public:
127 <  GeneralAtom(int theIndex): Atom(theIndex){}
127 >  GeneralAtom(int theIndex, SimState* theConfig): Atom(theIndex, theConfig){}
128    virtual ~GeneralAtom(){}
129  
130    int isDirectional( void ){ return 0; }
131 <  void zeroForces() {
218 <    frc[offsetX] = 0.0;
219 <    frc[offsetY] = 0.0;
220 <    frc[offsetZ] = 0.0;
221 <  }
131 >  void zeroForces( void );
132   };
133  
134   class DirectionalAtom : public Atom {
135    
136   public:
137 <  DirectionalAtom(int theIndex) : Atom(theIndex)
137 >  DirectionalAtom(int theIndex, SimState* theConfig) : Atom(theIndex,
138 >                                                            theConfig)
139    {
140      ssdIdentity = 0;
141      sux = 0.0;
142      suy = 0.0;
143      suz = 0.0;
144 +    myMu = 0.0;
145    }
146    virtual ~DirectionalAtom() {}
147  
148 +  virtual void setCoords(void);
149 +
150 +  void printAmatIndex( void );
151 +
152    int isDirectional(void) { return 1; }
153    
154    void setSSD( int value) { ssdIdentity = value; }
155    int isSSD(void) {return ssdIdentity; }
156  
241  void setA( double the_A[3][3] );
242
243  void setI( double the_I[3][3] );
244
245  void setQ( double the_q[4] );
157    
158    void setEuler( double phi, double theta, double psi );
159 +
160 +  double getSUx( void ) { return sux; }
161 +  double getSUy( void ) { return suy; }
162 +  double getSUz( void ) { return suz; }
163    
164    void setSUx( double the_sux ) { sux = the_sux; }
165    void setSUy( double the_suy ) { suy = the_suy; }
166    void setSUz( double the_suz ) { suz = the_suz; }
167  
168 <  void setJx( double the_jx ) { jx = the_jx; }
254 <  void setJy( double the_jy ) { jy = the_jy; }
255 <  void setJz( double the_jz ) { jz = the_jz; }
256 <    
257 <  void addTx( double the_tx ) { trq[offsetX] += the_tx;}
258 <  void addTy( double the_ty ) { trq[offsetY] += the_ty;}
259 <  void addTz( double the_tz ) { trq[offsetZ] += the_tz;}
168 >  void zeroForces();
169  
261  void zeroForces() {
262    frc[offsetX] = 0.0;
263    frc[offsetY] = 0.0;
264    frc[offsetZ] = 0.0;
265
266    trq[offsetX] = 0.0;
267    trq[offsetY] = 0.0;
268    trq[offsetZ] = 0.0;
269  }
270
271  double getAxx( void ) { return Amat[Axx]; }
272  double getAxy( void ) { return Amat[Axy]; }
273  double getAxz( void ) { return Amat[Axz]; }
274  
275  double getAyx( void ) { return Amat[Ayx]; }
276  double getAyy( void ) { return Amat[Ayy]; }
277  double getAyz( void ) { return Amat[Ayz]; }
278  
279  double getAzx( void ) { return Amat[Azx]; }
280  double getAzy( void ) { return Amat[Azy]; }
281  double getAzz( void ) { return Amat[Azz]; }
282
170    void getA( double the_A[3][3] ); // get the full rotation matrix
171 +  void setA( double the_A[3][3] );
172  
285  double getSUx( void ) { return sux; }
286  double getSUy( void ) { return suy; }
287  double getSUz( void ) { return suz; }
288
173    void getU( double the_u[3] ); // get the unit vetor
174 +  void updateU( void );
175 +
176    void getQ( double the_q[4] ); // get the quanternions
177 +  void setQ( double the_q[4] );
178  
179 +  void getJ( double theJ[3] );
180 +  void setJ( double theJ[3] );
181 +
182    double getJx( void ) { return jx; }
183    double getJy( void ) { return jy; }
184    double getJz( void ) { return jz; }
185  
186 <  double getTx( void ) { return trq[offsetX];}
187 <  double getTy( void ) { return trq[offsetY]; }
188 <  double getTz( void ) { return trq[offsetZ]; }
186 >  void setJx( double the_jx ) { jx = the_jx; }
187 >  void setJy( double the_jy ) { jy = the_jy; }
188 >  void setJz( double the_jz ) { jz = the_jz; }
189  
190 +  void getTrq( double theT[3] );
191 +  void addTrq( double theT[3] );
192 +
193 +  //  double getTx( void ) { return trq[offsetX];}
194 +  //  double getTy( void ) { return trq[offsetY]; }
195 +  //  double getTz( void ) { return trq[offsetZ]; }
196 +
197 +  void setI( double the_I[3][3] );
198 +  void getI( double the_I[3][3] );
199 +  
200    double getIxx( void ) { return Ixx; }
201    double getIxy( void ) { return Ixy; }
202    double getIxz( void ) { return Ixz; }
# Line 309 | Line 209 | class DirectionalAtom : public Atom { (public)
209    double getIzy( void ) { return Izy; }
210    double getIzz( void ) { return Izz; }
211  
212 <  double getMu( void ) { return mu[index]; }
213 <  void setMu( double the_mu ) { mu[index] = the_mu; }
212 >  double getMu( void );
213 >  void setMu( double the_mu );
214  
215    void lab2Body( double r[3] );
216    void body2Lab( double r[3] );
317  void updateU( void );
217  
218 +
219 +  // Four functions added for derivatives with respect to Euler Angles:
220 +  // (Needed for minimization routines):
221 +
222 +  void getGrad(double gradient[6] );
223 +  void getEulerAngles( double myEuler[3] );
224 +
225 +  double max(double x, double y);
226 +  double min(double x, double y);
227 +  
228 +
229   private:
230    int dIndex;
231  
232 +  double myMu;
233 +
234    double sux, suy, suz; // the standard unit vector    ( body fixed )
235    double jx, jy, jz;    // the angular momentum vector ( body fixed )
236    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines