mdtools/md_code/AllLong.cpp

#include <iostream>
#include <cstdlib>

#include "LRI.hpp"
#include "simError.h"

extern "C" {
  
  
  void long_range_force_(unsigned short int &exclude_check,
                         int* pair_i, int* pair_j, int &n_pairs, 
                         int &n_atoms, double *sigma, double *epslon, 
                         double &r_cut, double &r_list, 
                         double &box_x, double &box_y, double &box_z, 
                         double *R_x, double *R_y, double *R_z, 
                         double *F_x, double *F_y, double *F_z, 
                         double &potential_E,
                         unsigned short int &update, int *point,
                         int *list, int &maxnab, 
                         double *R_x0, double *R_y0, double *R_z0,
                         double* ux, double* uy, double* uz,
                         double* Tx, double* Ty, double* Tz,
                         double* Ex, double* Ey, double* Ez,
                         double* mu, double &rRF, double &rTaper, 
                         double &dielectric,
                         unsigned short int* isDipole,
                         unsigned short int* isVDW,
                         unsigned short int* isLJ,
                         unsigned short int* isSSD,
                         double* Axx, double *Axy, double *Axz,
                         double* Ayx, double *Ayy, double *Ayz,
                         double* Azx, double *Azy, double *Azz
                         );
  
  
  void long_range_check_(double &r_cut, double &r_list, 
                         unsigned short int &update,
                         int &n_atoms, double *R_x, double *R_y, double *R_z,
                         double *R_x0, double *R_y0, double *R_z0);
}


AllLong::AllLong( SimInfo* entry_plug ){
  int i, error;

  // get what info we need from the entry plug

  nPairs     = entry_plug->n_exclude;
  nAtoms     = entry_plug->n_atoms;
  atoms      = entry_plug->atoms;
  boxX       = entry_plug->box_x;
  boxY       = entry_plug->box_y;
  boxZ       = entry_plug->box_z;
  nDipoles   = entry_plug->n_dipoles;
  ex_pair *pair_list = entry_plug->excludes; 
  
  if( nDipoles ){
    dielectric = entry_plug->dielectric;
    rRF        = entry_plug->rRF;
    rTaper     = 0.95 * rRF;
  }
  else{
    dielectric = 0.0;
    rRF        = 0.0;
    rTaper     = 0.0;
  }
  
  exclude = 0;
  pairI = NULL;
  pairJ = NULL;

  if ( pair_list != NULL ){

    exclude = 1;
    
    pairI = new int[nPairs];
    pairJ = new int[nPairs];

    for(int i = 0; i < nPairs; i++){
    
      pairI[i] = pair_list[i].i;
      pairJ[i] = pair_list[i].j;
    }
  }

  sigma  = new double[nAtoms];
  epslon = new double[nAtoms];
  mu = new double[nAtoms];
  

  Rx  = new double[nAtoms];
  Ry  = new double[nAtoms];
  Rz  = new double[nAtoms];

  Rx0 = new double[nAtoms];
  Ry0 = new double[nAtoms];
  Rz0 = new double[nAtoms];

  Fx = new double[nAtoms];
  Fy = new double[nAtoms];
  Fz = new double[nAtoms];

  ux = new double[nAtoms];
  uy = new double[nAtoms];
  uz = new double[nAtoms];

  Tx = new double[nAtoms];
  Ty = new double[nAtoms];
  Tz = new double[nAtoms];

  Ex = new double[nAtoms];
  Ey = new double[nAtoms];
  Ez = new double[nAtoms];

  Axx = new double[nAtoms];
  Axy = new double[nAtoms];
  Axz = new double[nAtoms];

  Ayx = new double[nAtoms];
  Ayy = new double[nAtoms];
  Ayz = new double[nAtoms];

  Azx = new double[nAtoms];
  Azy = new double[nAtoms];
  Azz = new double[nAtoms];
  
  isDipole = new unsigned short int[nAtoms];
  isVDW    = new unsigned short int[nAtoms];
  isLJ     = new unsigned short int[nAtoms];
  isSSD    = new unsigned short int[nAtoms];

  // collect the array of dipole moments, mass, sigma, epsilon
  // also set the boolean arrays

  DirectionalAtom* dAtom;
  for( i=0; i<nAtoms; i++ ){
    
    isDipole[i] = 0;
    isVDW[i]    = 0;
    isLJ[i]     = 0;
    isSSD[i]    = 0;

    if( atoms[i]->hasDipole() ){

      dAtom = ( DirectionalAtom* )atoms[i];
      mu[i] = dAtom->getMu();
      isDipole[i] = 1;
      
      if( dAtom->isSSD() ) isSSD[i] = 1;
    }
    else mu[i] = 0.0;
    
    sigma[i]  = atoms[i]->getSigma();
    epslon[i] = atoms[i]->getEpslon();
    
    isVDW[i] = atoms[i]->isVDW();
    isLJ[i]  = atoms[i]->isLJ();
    
    // leech off the loop, and initialize the neighbor list initial positions

    Rx0[i] = 0.0; // these positions should force the first call to regenerate
    Ry0[i] = 0.0; // the neighbor list.
    Rz0[i] = 0.0;

    // leech again, and initialize the rotation matrixes

    Axx[i] = 1.0;
    Axy[i] = 0.0;
    Axz[i] = 0.0;

    Ayx[i] = 0.0;
    Ayy[i] = 1.0;
    Ayz[i] = 0.0;
    
    Azx[i] = 0.0;
    Azy[i] = 0.0;
    Azz[i] = 1.0;    
  }

  // check for force consistency

  error = 0;
  for( i=0; i<(nAtoms-1); i++ ){
    if( isVDW[i] && isLJ[i+1] ) error = 1;
    if( isLJ[i] && isVDW[i+1] ) error = 1;

    if( error ){
      sprintf(painCave.errMsg,
              "longRange force error: You cannot mix Lenard-Jones and"
              " \"Exp - r^6\" type interactions\n");
      painCave.isFatal = 1;
      simError();
    }
  }

  /* a little routine to set the rCut, should not exceed half of any
     box length */

  rCut = 0.0;
  double test_cut;
  
  for(int i = 0; i < nAtoms - 1; i++){
    for(int j = i+1; j < nAtoms; j++){
      
      test_cut = 2.5 * ( sigma[i] + sigma[j] ) / 2.0;
      if( test_cut > rCut) rCut = test_cut;
    }
  }

  if(rCut > (boxX / 2.0)) rCut = boxX / 2.0;
  if(rCut > (boxY / 2.0)) rCut = boxY / 2.0;
  if(rCut > (boxZ / 2.0)) rCut = boxZ / 2.0;

  if( rCut > rRF ){
    rList = rCut + 1.0;
    rListSmall = rCut;
  }
  else{
    rList = rRF + 1.0;
    rListSmall = rRF;
  }
  
  // the first time through, the neighbor list needs to be updated

  listUpdate = 1;

  // set the neighbor list info

  maxNab = nAtoms * 1000;
  point  = new int[nAtoms];
  list   = new int[maxNab];

  // give some love back to the entry plug

  if( entry_plug->longRange != NULL ) delete entry_plug->longRange;
  entry_plug->longRange = this;
}

AllLong::~AllLong(){

  if( pairI != NULL ) delete[] pairI;
  if( pairJ != NULL ) delete[] pairJ;

  delete[] sigma;
  delete[] epslon;
  delete[] mu;
  
  delete[] Rx;
  delete[] Ry;
  delete[] Rz;

  delete[] Rx0;
  delete[] Ry0;
  delete[] Rz0;

  delete[] Fx;
  delete[] Fy;
  delete[] Fz;

  delete[] ux;
  delete[] uy;
  delete[] uz;

  delete[] Tx;
  delete[] Ty;
  delete[] Tz;
  
  delete[] Ex;
  delete[] Ey;
  delete[] Ez; 

  delete[] Axx;
  delete[] Axy;
  delete[] Axz; 

  delete[] Ayx;
  delete[] Ayy;
  delete[] Ayz; 

  delete[] Azx;
  delete[] Azy;
  delete[] Azz; 
 
  delete[] isDipole;
  delete[] isVDW;
  delete[] isLJ;
  delete[] isSSD;

  delete[] point;
  delete[] list;
}


void AllLong::calc_forces( void ){
  
  int i;
  double u_temp[3];
  DirectionalAtom* dAtom;

  // fill our arrays;

  for( i=0; i<nAtoms; i++ ){

    if( isDipole[i] ){

      if( !(atoms[i]->isDirectional()) ){
        sprintf( painCave.errMsg,
                 "Something went Horribly wrong!!!!\n"
                 "We're all gonna die in here!!!!!!!\n"
                 "Oh yeah, by the way, atom %d" 
                 " somehow thinks it has a dipole,\n"
                 "but it isn't a directional atom type.\n", i);
        painCave.isFatal = 1;
        simError();
      }
               
      dAtom = ( DirectionalAtom* )atoms[i];
      dAtom->getU( u_temp );

      ux[i] = u_temp[0];
      uy[i] = u_temp[1];
      uz[i] = u_temp[2];

      Axx[i] = dAtom->getAxx();
      Axy[i] = dAtom->getAxy();
      Axz[i] = dAtom->getAxz();

      Ayx[i] = dAtom->getAyx();
      Ayy[i] = dAtom->getAyy();
      Ayz[i] = dAtom->getAyz();
      
      Azx[i] = dAtom->getAzx();
      Azy[i] = dAtom->getAzy();
      Azz[i] = dAtom->getAzz();

    }
    
    Rx[i] = atoms[i]->getX();
    Ry[i] = atoms[i]->getY();
    Rz[i] = atoms[i]->getZ();
  }


  // call fortran for the force calcs

  long_range_check_(rListSmall, rList, listUpdate,
                    nAtoms, Rx, Ry, Rz,
                    Rx0, Ry0, Rz0);

  long_range_force_(exclude,
                    pairI, pairJ, nPairs,
                    nAtoms, sigma, epslon, 
                    rCut, rList,
                    boxX, boxY, boxZ, 
                    Rx, Ry, Rz,
                    Fx, Fy, Fz, 
                    potentialE, 
                    listUpdate, point,
                    list, maxNab,
                    Rx0, Ry0, Rz0,
                    ux, uy, uz,
                    Tx, Ty, Tz,
                    Ex, Ey, Ez,
                    mu, rRF, rTaper,
                    dielectric,
                    isDipole, isVDW, isLJ, isSSD,
                    Axx, Axy, Axz,
                    Ayx, Ayy, Ayz,
                    Azx, Azy, Azz );

  // pass info back to the atoms

  for( i=0; i<nAtoms; i++ ){
    
    if( isDipole[i] ){
      
      dAtom = ( DirectionalAtom* )atoms[i];
      
      dAtom->addTx( Tx[i] );
      dAtom->addTy( Ty[i] );
      dAtom->addTz( Tz[i] );
    }

    atoms[i]->addFx( Fx[i] );
    atoms[i]->addFy( Fy[i] );
    atoms[i]->addFz( Fz[i] );
  }
}
Revision:	162
Committed:	Thu Oct 31 21:20:49 2002 UTC (21 years, 8 months ago) by mmeineke
File size:	8235 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	mmeineke	10	#include <iostream>
2			#include <cstdlib>
3
4			#include "LRI.hpp"
5	mmeineke	162	#include "simError.h"
6	mmeineke	10
7			extern "C" {
8
9
10			void long_range_force_(unsigned short int &exclude_check,
11			int* pair_i, int* pair_j, int &n_pairs,
12			int &n_atoms, double sigma, double epslon,
13			double &r_cut, double &r_list,
14			double &box_x, double &box_y, double &box_z,
15			double R_x, double R_y, double *R_z,
16			double F_x, double F_y, double *F_z,
17			double &potential_E,
18			unsigned short int &update, int *point,
19			int *list, int &maxnab,
20			double R_x0, double R_y0, double *R_z0,
21			double* ux, double* uy, double* uz,
22			double* Tx, double* Ty, double* Tz,
23			double* Ex, double* Ey, double* Ez,
24			double* mu, double &rRF, double &rTaper,
25			double &dielectric,
26			unsigned short int* isDipole,
27			unsigned short int* isVDW,
28			unsigned short int* isLJ,
29			unsigned short int* isSSD,
30			double* Axx, double Axy, double Axz,
31			double* Ayx, double Ayy, double Ayz,
32			double* Azx, double Azy, double Azz
33			);
34
35
36			void long_range_check_(double &r_cut, double &r_list,
37			unsigned short int &update,
38			int &n_atoms, double R_x, double R_y, double *R_z,
39			double R_x0, double R_y0, double *R_z0);
40			}
41
42
43
44			AllLong::AllLong( SimInfo* entry_plug ){
45			int i, error;
46
47			// get what info we need from the entry plug
48
49			nPairs = entry_plug->n_exclude;
50			nAtoms = entry_plug->n_atoms;
51			atoms = entry_plug->atoms;
52			boxX = entry_plug->box_x;
53			boxY = entry_plug->box_y;
54			boxZ = entry_plug->box_z;
55			nDipoles = entry_plug->n_dipoles;
56			ex_pair *pair_list = entry_plug->excludes;
57
58			if( nDipoles ){
59			dielectric = entry_plug->dielectric;
60			rRF = entry_plug->rRF;
61			rTaper = 0.95 * rRF;
62			}
63			else{
64			dielectric = 0.0;
65			rRF = 0.0;
66			rTaper = 0.0;
67			}
68
69			exclude = 0;
70			pairI = NULL;
71			pairJ = NULL;
72
73			if ( pair_list != NULL ){
74
75			exclude = 1;
76
77			pairI = new int[nPairs];
78			pairJ = new int[nPairs];
79
80			for(int i = 0; i < nPairs; i++){
81
82			pairI[i] = pair_list[i].i;
83			pairJ[i] = pair_list[i].j;
84			}
85			}
86
87			sigma = new double[nAtoms];
88			epslon = new double[nAtoms];
89			mu = new double[nAtoms];
90
91
92			Rx = new double[nAtoms];
93			Ry = new double[nAtoms];
94			Rz = new double[nAtoms];
95
96			Rx0 = new double[nAtoms];
97			Ry0 = new double[nAtoms];
98			Rz0 = new double[nAtoms];
99
100			Fx = new double[nAtoms];
101			Fy = new double[nAtoms];
102			Fz = new double[nAtoms];
103
104			ux = new double[nAtoms];
105			uy = new double[nAtoms];
106			uz = new double[nAtoms];
107
108			Tx = new double[nAtoms];
109			Ty = new double[nAtoms];
110			Tz = new double[nAtoms];
111
112			Ex = new double[nAtoms];
113			Ey = new double[nAtoms];
114			Ez = new double[nAtoms];
115
116			Axx = new double[nAtoms];
117			Axy = new double[nAtoms];
118			Axz = new double[nAtoms];
119
120			Ayx = new double[nAtoms];
121			Ayy = new double[nAtoms];
122			Ayz = new double[nAtoms];
123
124			Azx = new double[nAtoms];
125			Azy = new double[nAtoms];
126			Azz = new double[nAtoms];
127
128			isDipole = new unsigned short int[nAtoms];
129			isVDW = new unsigned short int[nAtoms];
130			isLJ = new unsigned short int[nAtoms];
131			isSSD = new unsigned short int[nAtoms];
132
133			// collect the array of dipole moments, mass, sigma, epsilon
134			// also set the boolean arrays
135
136			DirectionalAtom* dAtom;
137			for( i=0; i<nAtoms; i++ ){
138
139			isDipole[i] = 0;
140			isVDW[i] = 0;
141			isLJ[i] = 0;
142			isSSD[i] = 0;
143
144			if( atoms[i]->hasDipole() ){
145
146			dAtom = ( DirectionalAtom* )atoms[i];
147			mu[i] = dAtom->getMu();
148			isDipole[i] = 1;
149
150			if( dAtom->isSSD() ) isSSD[i] = 1;
151			}
152			else mu[i] = 0.0;
153
154			sigma[i] = atoms[i]->getSigma();
155			epslon[i] = atoms[i]->getEpslon();
156
157			isVDW[i] = atoms[i]->isVDW();
158			isLJ[i] = atoms[i]->isLJ();
159
160			// leech off the loop, and initialize the neighbor list initial positions
161
162			Rx0[i] = 0.0; // these positions should force the first call to regenerate
163			Ry0[i] = 0.0; // the neighbor list.
164			Rz0[i] = 0.0;
165
166			// leech again, and initialize the rotation matrixes
167
168			Axx[i] = 1.0;
169			Axy[i] = 0.0;
170			Axz[i] = 0.0;
171
172			Ayx[i] = 0.0;
173			Ayy[i] = 1.0;
174			Ayz[i] = 0.0;
175
176			Azx[i] = 0.0;
177			Azy[i] = 0.0;
178			Azz[i] = 1.0;
179			}
180
181			// check for force consistency
182
183			error = 0;
184			for( i=0; i<(nAtoms-1); i++ ){
185			if( isVDW[i] && isLJ[i+1] ) error = 1;
186			if( isLJ[i] && isVDW[i+1] ) error = 1;
187
188			if( error ){
189	mmeineke	162	sprintf(painCave.errMsg,
190			"longRange force error: You cannot mix Lenard-Jones and"
191			" \"Exp - r^6\" type interactions\n");
192			painCave.isFatal = 1;
193			simError();
194	mmeineke	10	}
195			}
196
197			/* a little routine to set the rCut, should not exceed half of any
198			box length */
199
200			rCut = 0.0;
201			double test_cut;
202
203			for(int i = 0; i < nAtoms - 1; i++){
204			for(int j = i+1; j < nAtoms; j++){
205
206			test_cut = 2.5 * ( sigma[i] + sigma[j] ) / 2.0;
207			if( test_cut > rCut) rCut = test_cut;
208			}
209			}
210
211			if(rCut > (boxX / 2.0)) rCut = boxX / 2.0;
212			if(rCut > (boxY / 2.0)) rCut = boxY / 2.0;
213			if(rCut > (boxZ / 2.0)) rCut = boxZ / 2.0;
214
215			if( rCut > rRF ){
216			rList = rCut + 1.0;
217			rListSmall = rCut;
218			}
219			else{
220			rList = rRF + 1.0;
221			rListSmall = rRF;
222			}
223
224			// the first time through, the neighbor list needs to be updated
225
226			listUpdate = 1;
227
228			// set the neighbor list info
229
230			maxNab = nAtoms * 1000;
231			point = new int[nAtoms];
232			list = new int[maxNab];
233
234			// give some love back to the entry plug
235
236			if( entry_plug->longRange != NULL ) delete entry_plug->longRange;
237			entry_plug->longRange = this;
238			}
239
240			AllLong::~AllLong(){
241
242			if( pairI != NULL ) delete[] pairI;
243			if( pairJ != NULL ) delete[] pairJ;
244
245			delete[] sigma;
246			delete[] epslon;
247			delete[] mu;
248
249			delete[] Rx;
250			delete[] Ry;
251			delete[] Rz;
252
253			delete[] Rx0;
254			delete[] Ry0;
255			delete[] Rz0;
256
257			delete[] Fx;
258			delete[] Fy;
259			delete[] Fz;
260
261			delete[] ux;
262			delete[] uy;
263			delete[] uz;
264
265			delete[] Tx;
266			delete[] Ty;
267			delete[] Tz;
268
269			delete[] Ex;
270			delete[] Ey;
271			delete[] Ez;
272
273			delete[] Axx;
274			delete[] Axy;
275			delete[] Axz;
276
277			delete[] Ayx;
278			delete[] Ayy;
279			delete[] Ayz;
280
281			delete[] Azx;
282			delete[] Azy;
283			delete[] Azz;
284
285			delete[] isDipole;
286			delete[] isVDW;
287			delete[] isLJ;
288			delete[] isSSD;
289
290			delete[] point;
291			delete[] list;
292			}
293
294
295			void AllLong::calc_forces( void ){
296
297			int i;
298			double u_temp[3];
299			DirectionalAtom* dAtom;
300
301			// fill our arrays;
302
303			for( i=0; i<nAtoms; i++ ){
304
305			if( isDipole[i] ){
306
307			if( !(atoms[i]->isDirectional()) ){
308	mmeineke	162	sprintf( painCave.errMsg,
309			"Something went Horribly wrong!!!!\n"
310			"We're all gonna die in here!!!!!!!\n"
311			"Oh yeah, by the way, atom %d"
312			" somehow thinks it has a dipole,\n"
313			"but it isn't a directional atom type.\n", i);
314			painCave.isFatal = 1;
315			simError();
316	mmeineke	10	}
317
318			dAtom = ( DirectionalAtom* )atoms[i];
319			dAtom->getU( u_temp );
320
321			ux[i] = u_temp[0];
322			uy[i] = u_temp[1];
323			uz[i] = u_temp[2];
324
325			Axx[i] = dAtom->getAxx();
326			Axy[i] = dAtom->getAxy();
327			Axz[i] = dAtom->getAxz();
328
329			Ayx[i] = dAtom->getAyx();
330			Ayy[i] = dAtom->getAyy();
331			Ayz[i] = dAtom->getAyz();
332
333			Azx[i] = dAtom->getAzx();
334			Azy[i] = dAtom->getAzy();
335			Azz[i] = dAtom->getAzz();
336
337			}
338
339			Rx[i] = atoms[i]->getX();
340			Ry[i] = atoms[i]->getY();
341			Rz[i] = atoms[i]->getZ();
342			}
343
344
345			// call fortran for the force calcs
346
347			long_range_check_(rListSmall, rList, listUpdate,
348			nAtoms, Rx, Ry, Rz,
349			Rx0, Ry0, Rz0);
350
351			long_range_force_(exclude,
352			pairI, pairJ, nPairs,
353			nAtoms, sigma, epslon,
354			rCut, rList,
355			boxX, boxY, boxZ,
356			Rx, Ry, Rz,
357			Fx, Fy, Fz,
358			potentialE,
359			listUpdate, point,
360			list, maxNab,
361			Rx0, Ry0, Rz0,
362			ux, uy, uz,
363			Tx, Ty, Tz,
364			Ex, Ey, Ez,
365			mu, rRF, rTaper,
366			dielectric,
367			isDipole, isVDW, isLJ, isSSD,
368			Axx, Axy, Axz,
369			Ayx, Ayy, Ayz,
370			Azx, Azy, Azz );
371
372			// pass info back to the atoms
373
374			for( i=0; i<nAtoms; i++ ){
375
376			if( isDipole[i] ){
377
378			dAtom = ( DirectionalAtom* )atoms[i];
379
380			dAtom->addTx( Tx[i] );
381			dAtom->addTy( Ty[i] );
382			dAtom->addTz( Tz[i] );
383			}
384
385			atoms[i]->addFx( Fx[i] );
386			atoms[i]->addFy( Fy[i] );
387			atoms[i]->addFz( Fz[i] );
388			}
389			}