ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/LJ_FF.cpp
Revision: 253
Committed: Thu Jan 30 15:20:21 2003 UTC (21 years, 5 months ago) by chuckv
File size: 14991 byte(s)
Log Message:
Added a generic util code directory and moved Linux_ifc_machdep to it.
MPI changes to compile MPI modules.

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