9 |
|
#include "SRI.hpp" |
10 |
|
#include "simError.h" |
11 |
|
|
12 |
– |
|
12 |
|
// Declare the structures that will be passed by the parser and MPI |
13 |
|
|
14 |
|
typedef struct{ |
24 |
|
int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ); |
25 |
|
|
26 |
|
#ifdef IS_MPI |
28 |
– |
|
27 |
|
#include "mpiForceField.h" |
28 |
|
|
29 |
|
MPI_Datatype mpiAtomStructType; |
31 |
|
#endif |
32 |
|
|
33 |
|
|
34 |
+ |
// declaration of functions needed to wrap the fortran module |
35 |
+ |
|
36 |
+ |
extern "C" { |
37 |
+ |
void forcefactory_( char* forceName, |
38 |
+ |
int* status, |
39 |
+ |
void (*wrapFunction)( void (*p1)( int* ident, |
40 |
+ |
double* mass, |
41 |
+ |
double* epslon, |
42 |
+ |
double* sigma, |
43 |
+ |
int* status ), |
44 |
+ |
void (*p2)( int *nLocal, |
45 |
+ |
int *identArray, |
46 |
+ |
int *isError ), |
47 |
+ |
void (*p3)( double* positionArray, |
48 |
+ |
double* forceArray, |
49 |
+ |
double* potentialEnergy, |
50 |
+ |
short int* doPotentialCalc )), |
51 |
+ |
int forceNameLength ); |
52 |
+ |
} |
53 |
+ |
|
54 |
+ |
|
55 |
+ |
void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, |
56 |
+ |
double* sigma, int* status ), |
57 |
+ |
void (*p2)( int *nLocal, int *identArray, int *isError ), |
58 |
+ |
void (*p3)( double* positionArray,double* forceArray, |
59 |
+ |
double* potentialEnergy, |
60 |
+ |
short int* doPotentialCalc ) ); |
61 |
+ |
|
62 |
+ |
void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma, |
63 |
+ |
int* status ); |
64 |
+ |
|
65 |
+ |
void (*initLJfortran) ( int *nLocal, int *identArray, int *isError ); |
66 |
+ |
|
67 |
+ |
LJ_FF* currentLJwrap; |
68 |
+ |
|
69 |
+ |
|
70 |
+ |
//**************************************************************** |
71 |
+ |
// begins the actual forcefield stuff. |
72 |
+ |
//**************************************************************** |
73 |
+ |
|
74 |
+ |
|
75 |
|
LJ_FF::LJ_FF(){ |
76 |
|
|
77 |
|
char fileName[200]; |
80 |
|
char temp[200]; |
81 |
|
char errMsg[1000]; |
82 |
|
|
83 |
+ |
// do the funtion wrapping |
84 |
+ |
currentLJwrap = this; |
85 |
+ |
wrapMe(); |
86 |
+ |
|
87 |
|
#ifdef IS_MPI |
88 |
|
int i; |
89 |
|
|
180 |
|
#endif // is_mpi |
181 |
|
} |
182 |
|
|
183 |
+ |
|
184 |
+ |
void LJ_FF::wrapMe( void ){ |
185 |
+ |
|
186 |
+ |
char* currentFF = "LJ"; |
187 |
+ |
int isError = 0; |
188 |
+ |
|
189 |
+ |
forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) ); |
190 |
+ |
|
191 |
+ |
if( isError ){ |
192 |
+ |
|
193 |
+ |
sprintf( painCave.errMsg, |
194 |
+ |
"LJ_FF error: an error was returned from fortran when the " |
195 |
+ |
"the functions were being wrapped.\n" ); |
196 |
+ |
painCave.isFatal = 1; |
197 |
+ |
simError(); |
198 |
+ |
} |
199 |
+ |
|
200 |
+ |
#ifdef IS_MPI |
201 |
+ |
sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." ); |
202 |
+ |
MPIcheckPoint(); |
203 |
+ |
#endif // is_mpi |
204 |
+ |
} |
205 |
+ |
|
206 |
+ |
|
207 |
+ |
void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, |
208 |
+ |
double* sigma, int* status ), |
209 |
+ |
void (*p2)( int*, int*, int* ), |
210 |
+ |
void (*p3)( double* positionArray,double* forceArray, |
211 |
+ |
double* potentialEnergy, |
212 |
+ |
short int* doPotentialCalc ) ){ |
213 |
+ |
|
214 |
+ |
|
215 |
+ |
newLJtype = p1; |
216 |
+ |
initLJfortran = p2; |
217 |
+ |
currentLJwrap->setLJfortran( p3 ); |
218 |
+ |
} |
219 |
+ |
|
220 |
+ |
|
221 |
+ |
|
222 |
|
void LJ_FF::initializeAtoms( void ){ |
223 |
|
|
224 |
|
class LinkedType { |
241 |
|
// check for duplicates |
242 |
|
|
243 |
|
if( !strcmp( info.name, name ) ){ |
244 |
< |
sprintf( simError.painCave, |
244 |
> |
sprintf( painCave.errMsg, |
245 |
|
"Duplicate LJ atom type \"%s\" found in " |
246 |
|
"the LJ_FF param file./n", |
247 |
|
name ); |
306 |
|
|
307 |
|
headAtomType = new LinkedType; |
308 |
|
|
309 |
< |
fastFoward( "AtomTypes", "initializeAtoms" ); |
309 |
> |
fastForward( "AtomTypes", "initializeAtoms" ); |
310 |
|
|
311 |
|
// we are now at the AtomTypes section. |
312 |
|
|
350 |
|
"LJ_FF atom structures read successfully." ); |
351 |
|
MPIcheckPoint(); |
352 |
|
|
353 |
< |
currentAtomType = headAtomType; |
353 |
> |
currentAtomType = headAtomType->next; //skip the first element who is a place holder. |
354 |
|
while( currentAtomType != NULL ){ |
355 |
|
currentAtomType->duplicate( info ); |
356 |
+ |
|
357 |
+ |
|
358 |
+ |
|
359 |
|
sendFrcStruct( &info, mpiAtomStructType ); |
360 |
+ |
|
361 |
+ |
sprintf( checkPointMsg, |
362 |
+ |
"successfully sent lJ force type: \"%s\"\n", |
363 |
+ |
info.name ); |
364 |
+ |
MPIcheckPoint(); |
365 |
+ |
|
366 |
|
currentAtomType = currentAtomType->next; |
367 |
|
} |
368 |
|
info.last = 1; |
378 |
|
|
379 |
|
headAtomType = new LinkedType; |
380 |
|
recieveFrcStruct( &info, mpiAtomStructType ); |
381 |
+ |
|
382 |
|
while( !info.last ){ |
383 |
|
|
384 |
+ |
|
385 |
+ |
|
386 |
|
headAtomType->add( info ); |
387 |
+ |
|
388 |
+ |
MPIcheckPoint(); |
389 |
+ |
|
390 |
|
recieveFrcStruct( &info, mpiAtomStructType ); |
391 |
|
} |
392 |
|
} |
393 |
|
#endif // is_mpi |
394 |
|
|
395 |
+ |
// call new A_types in fortran |
396 |
|
|
397 |
+ |
int isError; |
398 |
+ |
currentAtomType = headAtomType; |
399 |
+ |
while( currentAtomType != NULL ){ |
400 |
+ |
|
401 |
+ |
if( currentAtomType->name[0] != '\0' ){ |
402 |
+ |
isError = 0; |
403 |
+ |
newLJtype( &(currentAtomType->ident), |
404 |
+ |
&(currentAtomType->mass), |
405 |
+ |
&(currentAtomType->epslon), |
406 |
+ |
&(currentAtomType->sigma), |
407 |
+ |
&isError ); |
408 |
+ |
if( isError ){ |
409 |
+ |
sprintf( painCave.errMsg, |
410 |
+ |
"Error initializing the \"%s\" atom type in fortran\n", |
411 |
+ |
currentAtomType->name ); |
412 |
+ |
painCave.isFatal = 1; |
413 |
+ |
simError(); |
414 |
+ |
} |
415 |
+ |
} |
416 |
+ |
currentAtomType = currentAtomType->next; |
417 |
+ |
} |
418 |
+ |
|
419 |
+ |
#ifdef IS_MPI |
420 |
+ |
sprintf( checkPointMsg, |
421 |
+ |
"LJ_FF atom structures successfully sent to fortran\n" ); |
422 |
+ |
MPIcheckPoint(); |
423 |
+ |
#endif // is_mpi |
424 |
+ |
|
425 |
+ |
|
426 |
+ |
|
427 |
|
// initialize the atoms |
428 |
|
|
429 |
+ |
double bigSigma = 0.0; |
430 |
|
Atom* thisAtom; |
431 |
|
|
432 |
|
for( i=0; i<nAtoms; i++ ){ |
445 |
|
the_atoms[i]->setSigma( currentAtomType->sigma ); |
446 |
|
the_atoms[i]->setIdent( currentAtomType->ident ); |
447 |
|
the_atoms[i]->setLJ(); |
448 |
+ |
|
449 |
+ |
if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma; |
450 |
|
} |
451 |
|
|
452 |
+ |
|
453 |
+ |
#ifdef IS_MPI |
454 |
+ |
double tempBig = bigSigma; |
455 |
+ |
MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX ); |
456 |
+ |
#endif //is_mpi |
457 |
|
|
458 |
+ |
//calc rCut and rList |
459 |
+ |
|
460 |
+ |
entry_plug->rCut = 2.5 * bigSigma; |
461 |
+ |
if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0; |
462 |
+ |
if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0; |
463 |
+ |
if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0; |
464 |
+ |
|
465 |
+ |
entry_plug->rList = entry_plug->rCut + 1.0; |
466 |
+ |
|
467 |
|
// clean up the memory |
468 |
|
|
469 |
|
delete headAtomType; |
473 |
|
MPIcheckPoint(); |
474 |
|
#endif // is_mpi |
475 |
|
|
476 |
+ |
initFortran(); |
477 |
+ |
entry_plug->refreshSim(); |
478 |
+ |
|
479 |
|
} |
480 |
|
|
481 |
|
void LJ_FF::initializeBonds( bond_pair* the_bonds ){ |
617 |
|
} |
618 |
|
else return 0; |
619 |
|
} |
620 |
+ |
|
621 |
+ |
|
622 |
+ |
void LJ_FF::doForces( int calcPot ){ |
623 |
+ |
|
624 |
+ |
int i; |
625 |
+ |
double* frc; |
626 |
+ |
double* pos; |
627 |
+ |
short int passedCalcPot = (short int)calcPot; |
628 |
+ |
|
629 |
+ |
// forces are zeroed here, before any are acumulated. |
630 |
+ |
// NOTE: do not rezero the forces in Fortran. |
631 |
+ |
|
632 |
+ |
for(i=0; i<entry_plug->n_atoms; i++){ |
633 |
+ |
entry_plug->atoms[i]->zeroForces(); |
634 |
+ |
} |
635 |
+ |
|
636 |
+ |
frc = Atom::getFrcArray(); |
637 |
+ |
pos = Atom::getPosArray(); |
638 |
+ |
|
639 |
+ |
// entry_plug->lrPot = -1; |
640 |
+ |
doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot ); |
641 |
+ |
|
642 |
+ |
|
643 |
+ |
// fprintf( stderr, |
644 |
+ |
// "lrPot = %lf\n", entry_plug->lrPot ); |
645 |
+ |
|
646 |
+ |
} |
647 |
+ |
|
648 |
+ |
void LJ_FF::initFortran( void ){ |
649 |
+ |
|
650 |
+ |
int nLocal = entry_plug->n_atoms; |
651 |
+ |
int *ident; |
652 |
+ |
int isError; |
653 |
+ |
int i; |
654 |
+ |
|
655 |
+ |
ident = new int[nLocal]; |
656 |
+ |
|
657 |
+ |
for(i=0; i<nLocal; i++){ |
658 |
+ |
ident[i] = entry_plug->atoms[i]->getIdent(); |
659 |
+ |
} |
660 |
+ |
|
661 |
+ |
isError = 0; |
662 |
+ |
initLJfortran( &nLocal, ident, &isError ); |
663 |
+ |
|
664 |
+ |
if(isError){ |
665 |
+ |
sprintf( painCave.errMsg, |
666 |
+ |
"LJ_FF error: There was an error initializing the component list in fortran.\n" ); |
667 |
+ |
painCave.isFatal = 1; |
668 |
+ |
simError(); |
669 |
+ |
} |
670 |
+ |
|
671 |
+ |
|
672 |
+ |
#ifdef IS_MPI |
673 |
+ |
sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" ); |
674 |
+ |
MPIcheckPoint(); |
675 |
+ |
#endif // is_mpi |
676 |
+ |
|
677 |
+ |
delete[] ident; |
678 |
+ |
|
679 |
+ |
} |
680 |
+ |
|