ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/LJ_FF.cpp
Revision: 240
Committed: Wed Jan 22 21:45:20 2003 UTC (21 years, 5 months ago) by chuckv
File size: 14156 byte(s)
Log Message:
Added init function to c++ force module.

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)( void ),
211 void (*p3)( double* positionArray,double* forceArray,
212 double* potentialEnergy,
213 short int* doPotentialCalc ) ){
214
215
216 p1 = newLJtype;
217 p2 = initLJfortran;
218 this->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( simError.painCave,
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 fastFoward( "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 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
435
436 // clean up the memory
437
438 delete headAtomType;
439
440 #ifdef IS_MPI
441 sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
442 MPIcheckPoint();
443 #endif // is_mpi
444
445 initFortran();
446
447 }
448
449 void LJ_FF::initializeBonds( bond_pair* the_bonds ){
450
451 if( entry_plug->n_bonds ){
452 sprintf( painCave.errMsg,
453 "LJ_FF does not support bonds.\n" );
454 painCave.isFatal = 1;
455 simError();
456 }
457 #ifdef IS_MPI
458 MPIcheckPoint();
459 #endif // is_mpi
460
461 }
462
463 void LJ_FF::initializeBends( bend_set* the_bends ){
464
465 if( entry_plug->n_bends ){
466 sprintf( painCave.errMsg,
467 "LJ_FF does not support bends.\n" );
468 painCave.isFatal = 1;
469 simError();
470 }
471 #ifdef IS_MPI
472 MPIcheckPoint();
473 #endif // is_mpi
474
475 }
476
477 void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
478
479 if( entry_plug->n_torsions ){
480 sprintf( painCave.errMsg,
481 "LJ_FF does not support torsions.\n" );
482 painCave.isFatal = 1;
483 simError();
484 }
485 #ifdef IS_MPI
486 MPIcheckPoint();
487 #endif // is_mpi
488
489 }
490
491
492 void LJ_FF::fastForward( char* stopText, char* searchOwner ){
493
494 int foundText = 0;
495 char* the_token;
496
497 rewind( frcFile );
498 lineNum = 0;
499
500 eof_test = fgets( readLine, sizeof(readLine), frcFile );
501 lineNum++;
502 if( eof_test == NULL ){
503 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
504 " file is empty.\n",
505 searchOwner );
506 painCave.isFatal = 1;
507 simError();
508 }
509
510
511 while( !foundText ){
512 while( eof_test != NULL && readLine[0] != '#' ){
513 eof_test = fgets( readLine, sizeof(readLine), frcFile );
514 lineNum++;
515 }
516 if( eof_test == NULL ){
517 sprintf( painCave.errMsg,
518 "Error fast forwarding force file for %s at "
519 "line %d: file ended unexpectedly.\n",
520 searchOwner,
521 lineNum );
522 painCave.isFatal = 1;
523 simError();
524 }
525
526 the_token = strtok( readLine, " ,;\t#\n" );
527 foundText = !strcmp( stopText, the_token );
528
529 if( !foundText ){
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 }
544 }
545
546
547
548 int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ){
549
550 char* the_token;
551
552 the_token = strtok( lineBuffer, " \n\t,;" );
553 if( the_token != NULL ){
554
555 strcpy( info.name, the_token );
556
557 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
558 sprintf( painCave.errMsg,
559 "Error parseing AtomTypes: line %d\n", lineNum );
560 painCave.isFatal = 1;
561 simError();
562 }
563
564 info.mass = atof( the_token );
565
566 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
567 sprintf( painCave.errMsg,
568 "Error parseing AtomTypes: line %d\n", lineNum );
569 painCave.isFatal = 1;
570 simError();
571 }
572
573 info.epslon = atof( 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.sigma = atof( the_token );
583
584 return 1;
585 }
586 else return 0;
587 }
588
589
590 void LJ_FF::doForces( void ){
591
592 int i;
593 double* frc;
594 double* pos;
595 double potE;
596 short int calcPot = 0;
597
598 // forces are zeroed here, before any are acumulated.
599 // NOTE: do not rezero the forces in Fortran.
600
601 for(i=0; i<entry_plug->n_atoms; i++){
602 entry_plug->atoms[i]->zeroForces();
603 }
604
605 frc = Atom::getFrcArray();
606 pos = Atom::getPosArray();
607
608 doLJfortran( pos, frc, potE, calcPot );
609 }
610
611 void LJ_FF::initFortran( void ){
612
613 int nLocal = entry_plug->n_atoms;
614 int *ident;
615 int isError;
616 int i;
617
618 ident = new int[nLocal];
619
620 for(i=0; i<nLocal; i++){
621 ident[i] = entryplug->atoms[i]->getIdent();
622 }
623
624 isError = 0;
625 initLJfortran( &nLocal, ident, &isError );
626
627 if(isError){
628 sprintf( painCave.errMsg,
629 "LJ_FF error: There was an error initializing the component list in fortran.\n" );
630 painCave.isFatal = 1;
631 simError();
632 }
633
634
635 #ifdef IS_MPI
636 sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
637 MPIcheckPoint();
638 #endif // is_mpi
639
640 delete[] ident;
641
642 }
643