ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
Revision: 291
Committed: Wed Mar 5 20:35:54 2003 UTC (21 years, 4 months ago) by mmeineke
File size: 15631 byte(s)
Log Message:
overhauled TraPPE_Ex forcfield in C++

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