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 *
#	Content
1	#include <iostream>
2	#include <cstdlib>
3
4	#include "LRI.hpp"
5	#include "simError.h"
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	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	}
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	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	}
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	}