ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
Revision: 290
Committed: Thu Feb 27 21:25:47 2003 UTC (21 years, 6 months ago) by chuckv
File size: 15611 byte(s)
Log Message:
Made changes to lj_FF.cpp to pass stress tensor.

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