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, 4 months ago) by chuckv
File size: 15611 byte(s)
Log Message:
Made changes to lj_FF.cpp to pass stress tensor.

File Contents

# User Rev Content
1 mmeineke 270 #include <cstdlib>
2     #include <cstdio>
3     #include <cstring>
4    
5     #include <iostream>
6     using namespace std;
7    
8 mmeineke 275 #ifdef IS_MPI
9     #include <mpi.h>
10     #include <mpi++.h>
11     #endif //is_mpi
12    
13 mmeineke 270 #include "ForceFields.hpp"
14     #include "SRI.hpp"
15     #include "simError.h"
16    
17 mmeineke 275
18    
19 mmeineke 270 // 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 chuckv 290 double* tau,
58     short int* doPotentialCalc,
59     int* isError)),
60 mmeineke 270 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 chuckv 290 double* potentialEnergy, double* tau,
69     short int* doPotentialCalc, int* isError ) );
70 mmeineke 270
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 chuckv 290 int i, isError;
634 mmeineke 270 double* frc;
635     double* pos;
636 chuckv 290 double* tau;
637 mmeineke 270 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 chuckv 290 tau = entry_plug->tau;
649 mmeineke 270
650 chuckv 290 isError = 0;
651     doLJfortran( pos, frc, &(entry_plug->lrPot), tau, &passedCalcPot, &isError );
652 mmeineke 270
653    
654 chuckv 290 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 mmeineke 270 }
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