ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/LJ_FF.cpp
Revision: 263
Committed: Tue Feb 4 20:15:48 2003 UTC (21 years, 5 months ago) by chuckv
File size: 15216 byte(s)
Log Message:
Changed to single version.
No change.

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->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;
369 sendFrcStruct( &info, mpiAtomStructType );
370
371 }
372
373 else{
374
375 // listen for node 0 to send out the force params
376
377 MPIcheckPoint();
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++ ){
433
434 currentAtomType = headAtomType->find( the_atoms[i]->getType() );
435 if( currentAtomType == NULL ){
436 sprintf( painCave.errMsg,
437 "AtomType error, %s not found in force file.\n",
438 the_atoms[i]->getType() );
439 painCave.isFatal = 1;
440 simError();
441 }
442
443 the_atoms[i]->setMass( currentAtomType->mass );
444 the_atoms[i]->setEpslon( currentAtomType->epslon );
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;
470
471 #ifdef IS_MPI
472 sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
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 ){
482
483 if( entry_plug->n_bonds ){
484 sprintf( painCave.errMsg,
485 "LJ_FF does not support bonds.\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::initializeBends( bend_set* the_bends ){
496
497 if( entry_plug->n_bends ){
498 sprintf( painCave.errMsg,
499 "LJ_FF does not support bends.\n" );
500 painCave.isFatal = 1;
501 simError();
502 }
503 #ifdef IS_MPI
504 MPIcheckPoint();
505 #endif // is_mpi
506
507 }
508
509 void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
510
511 if( entry_plug->n_torsions ){
512 sprintf( painCave.errMsg,
513 "LJ_FF does not support torsions.\n" );
514 painCave.isFatal = 1;
515 simError();
516 }
517 #ifdef IS_MPI
518 MPIcheckPoint();
519 #endif // is_mpi
520
521 }
522
523
524 void LJ_FF::fastForward( char* stopText, char* searchOwner ){
525
526 int foundText = 0;
527 char* the_token;
528
529 rewind( frcFile );
530 lineNum = 0;
531
532 eof_test = fgets( readLine, sizeof(readLine), frcFile );
533 lineNum++;
534 if( eof_test == NULL ){
535 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
536 " file is empty.\n",
537 searchOwner );
538 painCave.isFatal = 1;
539 simError();
540 }
541
542
543 while( !foundText ){
544 while( eof_test != NULL && readLine[0] != '#' ){
545 eof_test = fgets( readLine, sizeof(readLine), frcFile );
546 lineNum++;
547 }
548 if( eof_test == NULL ){
549 sprintf( painCave.errMsg,
550 "Error fast forwarding force file for %s at "
551 "line %d: file ended unexpectedly.\n",
552 searchOwner,
553 lineNum );
554 painCave.isFatal = 1;
555 simError();
556 }
557
558 the_token = strtok( readLine, " ,;\t#\n" );
559 foundText = !strcmp( stopText, the_token );
560
561 if( !foundText ){
562 eof_test = fgets( readLine, sizeof(readLine), frcFile );
563 lineNum++;
564
565 if( eof_test == NULL ){
566 sprintf( painCave.errMsg,
567 "Error fast forwarding force file for %s at "
568 "line %d: file ended unexpectedly.\n",
569 searchOwner,
570 lineNum );
571 painCave.isFatal = 1;
572 simError();
573 }
574 }
575 }
576 }
577
578
579
580 int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ){
581
582 char* the_token;
583
584 the_token = strtok( lineBuffer, " \n\t,;" );
585 if( the_token != NULL ){
586
587 strcpy( info.name, the_token );
588
589 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
590 sprintf( painCave.errMsg,
591 "Error parseing AtomTypes: line %d\n", lineNum );
592 painCave.isFatal = 1;
593 simError();
594 }
595
596 info.mass = atof( the_token );
597
598 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
599 sprintf( painCave.errMsg,
600 "Error parseing AtomTypes: line %d\n", lineNum );
601 painCave.isFatal = 1;
602 simError();
603 }
604
605 info.epslon = atof( the_token );
606
607 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
608 sprintf( painCave.errMsg,
609 "Error parseing AtomTypes: line %d\n", lineNum );
610 painCave.isFatal = 1;
611 simError();
612 }
613
614 info.sigma = atof( the_token );
615
616 return 1;
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