ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/LJ_FF.cpp
Revision: 254
Committed: Thu Jan 30 20:03:37 2003 UTC (21 years, 5 months ago) by chuckv
File size: 14989 byte(s)
Log Message:
Bug fixes for mpi version of code

File Contents

# Content
1 #include <cstdlib>
2 #include <cstdio>
3 #include <cstring>
4
5 #include <iostream>
6 using namespace std;
7
8 #include "ForceFields.hpp"
9 #include "SRI.hpp"
10 #include "simError.h"
11
12 // Declare the structures that will be passed by the parser and MPI
13
14 typedef struct{
15 char name[15];
16 double mass;
17 double epslon;
18 double sigma;
19 int ident;
20 int last; // 0 -> default
21 // 1 -> in MPI: tells nodes to stop listening
22 } atomStruct;
23
24 int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info );
25
26 #ifdef IS_MPI
27 #include "mpiForceField.h"
28
29 MPI_Datatype mpiAtomStructType;
30
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];
78 char* ffPath_env = "FORCE_PARAM_PATH";
79 char* ffPath;
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
90 // **********************************************************************
91 // Init the atomStruct mpi type
92
93 atomStruct atomProto; // mpiPrototype
94 int atomBC[3] = {15,3,2}; // block counts
95 MPI_Aint atomDspls[3]; // displacements
96 MPI_Datatype atomMbrTypes[3]; // member mpi types
97
98 MPI_Address(&atomProto.name, &atomDspls[0]);
99 MPI_Address(&atomProto.mass, &atomDspls[1]);
100 MPI_Address(&atomProto.ident, &atomDspls[2]);
101
102 atomMbrTypes[0] = MPI_CHAR;
103 atomMbrTypes[1] = MPI_DOUBLE;
104 atomMbrTypes[2] = MPI_INT;
105
106 for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
107
108 MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
109 MPI_Type_commit(&mpiAtomStructType);
110
111 // ***********************************************************************
112
113 if( worldRank == 0 ){
114 #endif
115
116 // generate the force file name
117
118 strcpy( fileName, "LJ_FF.frc" );
119 // fprintf( stderr,"Trying to open %s\n", fileName );
120
121 // attempt to open the file in the current directory first.
122
123 frcFile = fopen( fileName, "r" );
124
125 if( frcFile == NULL ){
126
127 // next see if the force path enviorment variable is set
128
129 ffPath = getenv( ffPath_env );
130 if( ffPath == NULL ) {
131 sprintf( painCave.errMsg,
132 "Error opening the force field parameter file: %s\n"
133 "Have you tried setting the FORCE_PARAM_PATH environment "
134 "vairable?\n",
135 fileName );
136 painCave.isFatal = 1;
137 simError();
138 }
139
140
141 strcpy( temp, ffPath );
142 strcat( temp, "/" );
143 strcat( temp, fileName );
144 strcpy( fileName, temp );
145
146 frcFile = fopen( fileName, "r" );
147
148 if( frcFile == NULL ){
149
150 sprintf( painCave.errMsg,
151 "Error opening the force field parameter file: %s\n"
152 "Have you tried setting the FORCE_PARAM_PATH environment "
153 "vairable?\n",
154 fileName );
155 painCave.isFatal = 1;
156 simError();
157 }
158 }
159
160 #ifdef IS_MPI
161 }
162
163 sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
164 MPIcheckPoint();
165
166 #endif // is_mpi
167 }
168
169
170 LJ_FF::~LJ_FF(){
171
172 #ifdef IS_MPI
173 if( worldRank == 0 ){
174 #endif // is_mpi
175
176 fclose( frcFile );
177
178 #ifdef IS_MPI
179 }
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 {
225 public:
226 LinkedType(){
227 next = NULL;
228 name[0] = '\0';
229 }
230 ~LinkedType(){ if( next != NULL ) delete next; }
231
232 LinkedType* find(char* key){
233 if( !strcmp(name, key) ) return this;
234 if( next != NULL ) return next->find(key);
235 return NULL;
236 }
237
238
239 void add( atomStruct &info ){
240
241 // check for duplicates
242
243 if( !strcmp( info.name, name ) ){
244 sprintf( painCave.errMsg,
245 "Duplicate LJ atom type \"%s\" found in "
246 "the LJ_FF param file./n",
247 name );
248 painCave.isFatal = 1;
249 simError();
250 }
251
252 if( next != NULL ) next->add(info);
253 else{
254 next = new LinkedType();
255 strcpy(next->name, info.name);
256 next->mass = info.mass;
257 next->epslon = info.epslon;
258 next->sigma = info.sigma;
259 next->ident = info.ident;
260 }
261 }
262
263
264 #ifdef IS_MPI
265
266 void duplicate( atomStruct &info ){
267 strcpy(info.name, name);
268 info.mass = mass;
269 info.epslon = epslon;
270 info.sigma = sigma;
271 info.ident = ident;
272 info.last = 0;
273 }
274
275
276 #endif
277
278 char name[15];
279 double mass;
280 double epslon;
281 double sigma;
282 int ident;
283 LinkedType* next;
284 };
285
286 LinkedType* headAtomType;
287 LinkedType* currentAtomType;
288 atomStruct info;
289 info.last = 1; // initialize last to have the last set.
290 // if things go well, last will be set to 0
291
292 int i;
293 int identNum;
294
295 Atom** the_atoms;
296 int nAtoms;
297 the_atoms = entry_plug->atoms;
298 nAtoms = entry_plug->n_atoms;
299
300
301 #ifdef IS_MPI
302 if( worldRank == 0 ){
303 #endif
304
305 // read in the atom types.
306
307 headAtomType = new LinkedType;
308
309 fastForward( "AtomTypes", "initializeAtoms" );
310
311 // we are now at the AtomTypes section.
312
313 eof_test = fgets( readLine, sizeof(readLine), frcFile );
314 lineNum++;
315
316
317 // read a line, and start parseing out the atom types
318
319 if( eof_test == NULL ){
320 sprintf( painCave.errMsg,
321 "Error in reading Atoms from force file at line %d.\n",
322 lineNum );
323 painCave.isFatal = 1;
324 simError();
325 }
326
327 identNum = 1;
328 // stop reading at end of file, or at next section
329 while( readLine[0] != '#' && eof_test != NULL ){
330
331 // toss comment lines
332 if( readLine[0] != '!' ){
333
334 // the parser returns 0 if the line was blank
335 if( parseAtomLJ( readLine, lineNum, info ) ){
336 info.ident = identNum;
337 headAtomType->add( info );;
338 identNum++;
339 }
340 }
341 eof_test = fgets( readLine, sizeof(readLine), frcFile );
342 lineNum++;
343 }
344
345 #ifdef IS_MPI
346
347 // send out the linked list to all the other processes
348
349 sprintf( checkPointMsg,
350 "LJ_FF atom structures read successfully." );
351 MPIcheckPoint();
352
353 currentAtomType = headAtomType;
354 while( currentAtomType != NULL ){
355 currentAtomType->duplicate( info );
356 sendFrcStruct( &info, mpiAtomStructType );
357 currentAtomType = currentAtomType->next;
358 }
359 info.last = 1;
360 sendFrcStruct( &info, mpiAtomStructType );
361
362 }
363
364 else{
365
366 // listen for node 0 to send out the force params
367
368 MPIcheckPoint();
369
370 headAtomType = new LinkedType;
371 recieveFrcStruct( &info, mpiAtomStructType );
372 while( !info.last ){
373
374 headAtomType->add( info );
375 recieveFrcStruct( &info, mpiAtomStructType );
376 }
377 }
378 #endif // is_mpi
379
380 // call new A_types in fortran
381
382 int isError;
383 currentAtomType = headAtomType;
384 while( currentAtomType != NULL ){
385
386 if( currentAtomType->name[0] != '\0' ){
387 isError = 0;
388 newLJtype( &(currentAtomType->ident),
389 &(currentAtomType->mass),
390 &(currentAtomType->epslon),
391 &(currentAtomType->sigma),
392 &isError );
393 if( isError ){
394 sprintf( painCave.errMsg,
395 "Error initializing the \"%s\" atom type in fortran\n",
396 currentAtomType->name );
397 painCave.isFatal = 1;
398 simError();
399 }
400 }
401 currentAtomType = currentAtomType->next;
402 }
403
404 #ifdef IS_MPI
405 sprintf( checkPointMsg,
406 "LJ_FF atom structures successfully sent to fortran\n" );
407 MPIcheckPoint();
408 #endif // is_mpi
409
410
411
412 // initialize the atoms
413
414 double bigSigma = 0.0;
415 Atom* thisAtom;
416
417 for( i=0; i<nAtoms; i++ ){
418
419 currentAtomType = headAtomType->find( the_atoms[i]->getType() );
420 if( currentAtomType == NULL ){
421 sprintf( painCave.errMsg,
422 "AtomType error, %s not found in force file.\n",
423 the_atoms[i]->getType() );
424 painCave.isFatal = 1;
425 simError();
426 }
427
428 the_atoms[i]->setMass( currentAtomType->mass );
429 the_atoms[i]->setEpslon( currentAtomType->epslon );
430 the_atoms[i]->setSigma( currentAtomType->sigma );
431 the_atoms[i]->setIdent( currentAtomType->ident );
432 the_atoms[i]->setLJ();
433
434 if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
435 }
436
437
438 #ifdef IS_MPI
439 double tempBig = bigSigma;
440 MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
441 #endif //is_mpi
442
443 //calc rCut and rList
444
445 entry_plug->rCut = 2.5 * bigSigma;
446 if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
447 if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
448 if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;
449
450 entry_plug->rList = entry_plug->rCut + 1.0;
451
452 // clean up the memory
453
454 delete headAtomType;
455
456 #ifdef IS_MPI
457 sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
458 MPIcheckPoint();
459 #endif // is_mpi
460
461 initFortran();
462 entry_plug->refreshSim();
463
464 }
465
466 void LJ_FF::initializeBonds( bond_pair* the_bonds ){
467
468 if( entry_plug->n_bonds ){
469 sprintf( painCave.errMsg,
470 "LJ_FF does not support bonds.\n" );
471 painCave.isFatal = 1;
472 simError();
473 }
474 #ifdef IS_MPI
475 MPIcheckPoint();
476 #endif // is_mpi
477
478 }
479
480 void LJ_FF::initializeBends( bend_set* the_bends ){
481
482 if( entry_plug->n_bends ){
483 sprintf( painCave.errMsg,
484 "LJ_FF does not support bends.\n" );
485 painCave.isFatal = 1;
486 simError();
487 }
488 #ifdef IS_MPI
489 MPIcheckPoint();
490 #endif // is_mpi
491
492 }
493
494 void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
495
496 if( entry_plug->n_torsions ){
497 sprintf( painCave.errMsg,
498 "LJ_FF does not support torsions.\n" );
499 painCave.isFatal = 1;
500 simError();
501 }
502 #ifdef IS_MPI
503 MPIcheckPoint();
504 #endif // is_mpi
505
506 }
507
508
509 void LJ_FF::fastForward( char* stopText, char* searchOwner ){
510
511 int foundText = 0;
512 char* the_token;
513
514 rewind( frcFile );
515 lineNum = 0;
516
517 eof_test = fgets( readLine, sizeof(readLine), frcFile );
518 lineNum++;
519 if( eof_test == NULL ){
520 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
521 " file is empty.\n",
522 searchOwner );
523 painCave.isFatal = 1;
524 simError();
525 }
526
527
528 while( !foundText ){
529 while( eof_test != NULL && readLine[0] != '#' ){
530 eof_test = fgets( readLine, sizeof(readLine), frcFile );
531 lineNum++;
532 }
533 if( eof_test == NULL ){
534 sprintf( painCave.errMsg,
535 "Error fast forwarding force file for %s at "
536 "line %d: file ended unexpectedly.\n",
537 searchOwner,
538 lineNum );
539 painCave.isFatal = 1;
540 simError();
541 }
542
543 the_token = strtok( readLine, " ,;\t#\n" );
544 foundText = !strcmp( stopText, the_token );
545
546 if( !foundText ){
547 eof_test = fgets( readLine, sizeof(readLine), frcFile );
548 lineNum++;
549
550 if( eof_test == NULL ){
551 sprintf( painCave.errMsg,
552 "Error fast forwarding force file for %s at "
553 "line %d: file ended unexpectedly.\n",
554 searchOwner,
555 lineNum );
556 painCave.isFatal = 1;
557 simError();
558 }
559 }
560 }
561 }
562
563
564
565 int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ){
566
567 char* the_token;
568
569 the_token = strtok( lineBuffer, " \n\t,;" );
570 if( the_token != NULL ){
571
572 strcpy( info.name, the_token );
573
574 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
575 sprintf( painCave.errMsg,
576 "Error parseing AtomTypes: line %d\n", lineNum );
577 painCave.isFatal = 1;
578 simError();
579 }
580
581 info.mass = atof( the_token );
582
583 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
584 sprintf( painCave.errMsg,
585 "Error parseing AtomTypes: line %d\n", lineNum );
586 painCave.isFatal = 1;
587 simError();
588 }
589
590 info.epslon = atof( the_token );
591
592 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
593 sprintf( painCave.errMsg,
594 "Error parseing AtomTypes: line %d\n", lineNum );
595 painCave.isFatal = 1;
596 simError();
597 }
598
599 info.sigma = atof( the_token );
600
601 return 1;
602 }
603 else return 0;
604 }
605
606
607 void LJ_FF::doForces( int calcPot ){
608
609 int i;
610 double* frc;
611 double* pos;
612 short int passedCalcPot = (short int)calcPot;
613
614 // forces are zeroed here, before any are acumulated.
615 // NOTE: do not rezero the forces in Fortran.
616
617 for(i=0; i<entry_plug->n_atoms; i++){
618 entry_plug->atoms[i]->zeroForces();
619 }
620
621 frc = Atom::getFrcArray();
622 pos = Atom::getPosArray();
623
624 // entry_plug->lrPot = -1;
625 doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot );
626
627
628 // fprintf( stderr,
629 // "lrPot = %lf\n", entry_plug->lrPot );
630
631 }
632
633 void LJ_FF::initFortran( void ){
634
635 int nLocal = entry_plug->n_atoms;
636 int *ident;
637 int isError;
638 int i;
639
640 ident = new int[nLocal];
641
642 for(i=0; i<nLocal; i++){
643 ident[i] = entry_plug->atoms[i]->getIdent();
644 }
645
646 isError = 0;
647 initLJfortran( &nLocal, ident, &isError );
648
649 if(isError){
650 sprintf( painCave.errMsg,
651 "LJ_FF error: There was an error initializing the component list in fortran.\n" );
652 painCave.isFatal = 1;
653 simError();
654 }
655
656
657 #ifdef IS_MPI
658 sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
659 MPIcheckPoint();
660 #endif // is_mpi
661
662 delete[] ident;
663
664 }
665