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:	189
Committed:	Tue Nov 26 21:04:43 2002 UTC (22 years, 11 months ago) by mmeineke
File size:	8236 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
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	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	}
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	}