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 exrListceed 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:	248
Committed:	Mon Jan 27 19:28:21 2003 UTC (21 years, 5 months ago) by chuckv
File size:	8241 byte(s)
Log Message:	final version before the single processor build
#	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	chuckv	248	/* a little routine to set the rCut, should not exrListceed half of any
198	mmeineke	10	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	mmeineke	189
239	mmeineke	10	}
240
241			AllLong::~AllLong(){
242
243			if( pairI != NULL ) delete[] pairI;
244			if( pairJ != NULL ) delete[] pairJ;
245
246			delete[] sigma;
247			delete[] epslon;
248			delete[] mu;
249
250			delete[] Rx;
251			delete[] Ry;
252			delete[] Rz;
253
254			delete[] Rx0;
255			delete[] Ry0;
256			delete[] Rz0;
257
258			delete[] Fx;
259			delete[] Fy;
260			delete[] Fz;
261
262			delete[] ux;
263			delete[] uy;
264			delete[] uz;
265
266			delete[] Tx;
267			delete[] Ty;
268			delete[] Tz;
269
270			delete[] Ex;
271			delete[] Ey;
272			delete[] Ez;
273
274			delete[] Axx;
275			delete[] Axy;
276			delete[] Axz;
277
278			delete[] Ayx;
279			delete[] Ayy;
280			delete[] Ayz;
281
282			delete[] Azx;
283			delete[] Azy;
284			delete[] Azz;
285
286			delete[] isDipole;
287			delete[] isVDW;
288			delete[] isLJ;
289			delete[] isSSD;
290
291			delete[] point;
292			delete[] list;
293			}
294
295
296			void AllLong::calc_forces( void ){
297
298			int i;
299			double u_temp[3];
300			DirectionalAtom* dAtom;
301
302			// fill our arrays;
303
304			for( i=0; i<nAtoms; i++ ){
305
306			if( isDipole[i] ){
307
308			if( !(atoms[i]->isDirectional()) ){
309	mmeineke	162	sprintf( painCave.errMsg,
310			"Something went Horribly wrong!!!!\n"
311			"We're all gonna die in here!!!!!!!\n"
312			"Oh yeah, by the way, atom %d"
313			" somehow thinks it has a dipole,\n"
314			"but it isn't a directional atom type.\n", i);
315			painCave.isFatal = 1;
316			simError();
317	mmeineke	10	}
318
319			dAtom = ( DirectionalAtom* )atoms[i];
320			dAtom->getU( u_temp );
321
322			ux[i] = u_temp[0];
323			uy[i] = u_temp[1];
324			uz[i] = u_temp[2];
325
326			Axx[i] = dAtom->getAxx();
327			Axy[i] = dAtom->getAxy();
328			Axz[i] = dAtom->getAxz();
329
330			Ayx[i] = dAtom->getAyx();
331			Ayy[i] = dAtom->getAyy();
332			Ayz[i] = dAtom->getAyz();
333
334			Azx[i] = dAtom->getAzx();
335			Azy[i] = dAtom->getAzy();
336			Azz[i] = dAtom->getAzz();
337
338			}
339
340			Rx[i] = atoms[i]->getX();
341			Ry[i] = atoms[i]->getY();
342			Rz[i] = atoms[i]->getZ();
343			}
344
345
346			// call fortran for the force calcs
347
348			long_range_check_(rListSmall, rList, listUpdate,
349			nAtoms, Rx, Ry, Rz,
350			Rx0, Ry0, Rz0);
351
352			long_range_force_(exclude,
353			pairI, pairJ, nPairs,
354			nAtoms, sigma, epslon,
355			rCut, rList,
356			boxX, boxY, boxZ,
357			Rx, Ry, Rz,
358			Fx, Fy, Fz,
359			potentialE,
360			listUpdate, point,
361			list, maxNab,
362			Rx0, Ry0, Rz0,
363			ux, uy, uz,
364			Tx, Ty, Tz,
365			Ex, Ey, Ez,
366			mu, rRF, rTaper,
367			dielectric,
368			isDipole, isVDW, isLJ, isSSD,
369			Axx, Axy, Axz,
370			Ayx, Ayy, Ayz,
371			Azx, Azy, Azz );
372
373			// pass info back to the atoms
374
375			for( i=0; i<nAtoms; i++ ){
376
377			if( isDipole[i] ){
378
379			dAtom = ( DirectionalAtom* )atoms[i];
380
381			dAtom->addTx( Tx[i] );
382			dAtom->addTy( Ty[i] );
383			dAtom->addTz( Tz[i] );
384			}
385
386			atoms[i]->addFx( Fx[i] );
387			atoms[i]->addFy( Fy[i] );
388			atoms[i]->addFz( Fz[i] );
389			}
390			}