ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
Revision: 275
Committed: Tue Feb 18 21:06:36 2003 UTC (21 years, 5 months ago) by mmeineke
File size: 15285 byte(s)
Log Message:
libmdCode builds.

File Contents

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