mdtools/md_code/AllLong.cpp

#include <iostream>
#include <cstdlib>

#include "LRI.hpp"


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 ){
      std::cerr << "longRange force error: You cannot mix Lenard-Jones and"
                << " \"Exp - r^6\" type interactions\n";
      exit(8);
    }
  }

  /* 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()) ){
        std::cerr << "Something went Horribly wrong!!!!\n"
                  << "We're all gonna die in here!!!!!!!\n"
                  << "Oh yeah, by the way, atom " << i 
                  << " somehow thinks it has a dipole,\n"
                  << "but it isn't a directional atom type.\n";
        exit(8);
      }
               
      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:	10
Committed:	Tue Jul 9 18:40:59 2002 UTC (22 years ago) by mmeineke
Original Path:	*branches/mmeineke/mdtools/md_code/AllLong.cpp*
File size:	8106 byte(s)
Log Message:	everything you need to make libmdtools
#	User	Rev	Content
1	mmeineke	10	#include <iostream>
2			#include <cstdlib>
3
4			#include "LRI.hpp"
5
6
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			std::cerr << "longRange force error: You cannot mix Lenard-Jones and"
190			<< " \"Exp - r^6\" type interactions\n";
191			exit(8);
192			}
193			}
194
195			/* a little routine to set the rCut, should not exceed half of any
196			box length */
197
198			rCut = 0.0;
199			double test_cut;
200
201			for(int i = 0; i < nAtoms - 1; i++){
202			for(int j = i+1; j < nAtoms; j++){
203
204			test_cut = 2.5 * ( sigma[i] + sigma[j] ) / 2.0;
205			if( test_cut > rCut) rCut = test_cut;
206			}
207			}
208
209			if(rCut > (boxX / 2.0)) rCut = boxX / 2.0;
210			if(rCut > (boxY / 2.0)) rCut = boxY / 2.0;
211			if(rCut > (boxZ / 2.0)) rCut = boxZ / 2.0;
212
213			if( rCut > rRF ){
214			rList = rCut + 1.0;
215			rListSmall = rCut;
216			}
217			else{
218			rList = rRF + 1.0;
219			rListSmall = rRF;
220			}
221
222			// the first time through, the neighbor list needs to be updated
223
224			listUpdate = 1;
225
226			// set the neighbor list info
227
228			maxNab = nAtoms * 1000;
229			point = new int[nAtoms];
230			list = new int[maxNab];
231
232			// give some love back to the entry plug
233
234			if( entry_plug->longRange != NULL ) delete entry_plug->longRange;
235			entry_plug->longRange = this;
236			}
237
238			AllLong::~AllLong(){
239
240			if( pairI != NULL ) delete[] pairI;
241			if( pairJ != NULL ) delete[] pairJ;
242
243			delete[] sigma;
244			delete[] epslon;
245			delete[] mu;
246
247			delete[] Rx;
248			delete[] Ry;
249			delete[] Rz;
250
251			delete[] Rx0;
252			delete[] Ry0;
253			delete[] Rz0;
254
255			delete[] Fx;
256			delete[] Fy;
257			delete[] Fz;
258
259			delete[] ux;
260			delete[] uy;
261			delete[] uz;
262
263			delete[] Tx;
264			delete[] Ty;
265			delete[] Tz;
266
267			delete[] Ex;
268			delete[] Ey;
269			delete[] Ez;
270
271			delete[] Axx;
272			delete[] Axy;
273			delete[] Axz;
274
275			delete[] Ayx;
276			delete[] Ayy;
277			delete[] Ayz;
278
279			delete[] Azx;
280			delete[] Azy;
281			delete[] Azz;
282
283			delete[] isDipole;
284			delete[] isVDW;
285			delete[] isLJ;
286			delete[] isSSD;
287
288			delete[] point;
289			delete[] list;
290			}
291
292
293			void AllLong::calc_forces( void ){
294
295			int i;
296			double u_temp[3];
297			DirectionalAtom* dAtom;
298
299			// fill our arrays;
300
301			for( i=0; i<nAtoms; i++ ){
302
303			if( isDipole[i] ){
304
305			if( !(atoms[i]->isDirectional()) ){
306			std::cerr << "Something went Horribly wrong!!!!\n"
307			<< "We're all gonna die in here!!!!!!!\n"
308			<< "Oh yeah, by the way, atom " << i
309			<< " somehow thinks it has a dipole,\n"
310			<< "but it isn't a directional atom type.\n";
311			exit(8);
312			}
313
314			dAtom = ( DirectionalAtom* )atoms[i];
315			dAtom->getU( u_temp );
316
317			ux[i] = u_temp[0];
318			uy[i] = u_temp[1];
319			uz[i] = u_temp[2];
320
321			Axx[i] = dAtom->getAxx();
322			Axy[i] = dAtom->getAxy();
323			Axz[i] = dAtom->getAxz();
324
325			Ayx[i] = dAtom->getAyx();
326			Ayy[i] = dAtom->getAyy();
327			Ayz[i] = dAtom->getAyz();
328
329			Azx[i] = dAtom->getAzx();
330			Azy[i] = dAtom->getAzy();
331			Azz[i] = dAtom->getAzz();
332
333			}
334
335			Rx[i] = atoms[i]->getX();
336			Ry[i] = atoms[i]->getY();
337			Rz[i] = atoms[i]->getZ();
338			}
339
340
341			// call fortran for the force calcs
342
343			long_range_check_(rListSmall, rList, listUpdate,
344			nAtoms, Rx, Ry, Rz,
345			Rx0, Ry0, Rz0);
346
347			long_range_force_(exclude,
348			pairI, pairJ, nPairs,
349			nAtoms, sigma, epslon,
350			rCut, rList,
351			boxX, boxY, boxZ,
352			Rx, Ry, Rz,
353			Fx, Fy, Fz,
354			potentialE,
355			listUpdate, point,
356			list, maxNab,
357			Rx0, Ry0, Rz0,
358			ux, uy, uz,
359			Tx, Ty, Tz,
360			Ex, Ey, Ez,
361			mu, rRF, rTaper,
362			dielectric,
363			isDipole, isVDW, isLJ, isSSD,
364			Axx, Axy, Axz,
365			Ayx, Ayy, Ayz,
366			Azx, Azy, Azz );
367
368			// pass info back to the atoms
369
370			for( i=0; i<nAtoms; i++ ){
371
372			if( isDipole[i] ){
373
374			dAtom = ( DirectionalAtom* )atoms[i];
375
376			dAtom->addTx( Tx[i] );
377			dAtom->addTy( Ty[i] );
378			dAtom->addTz( Tz[i] );
379			}
380
381			atoms[i]->addFx( Fx[i] );
382			atoms[i]->addFy( Fy[i] );
383			atoms[i]->addFz( Fz[i] );
384			}
385			}