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
File size:	8106 byte(s)
Log Message:	everything you need to make libmdtools
#	Content
1	#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	}