ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
Revision: 321
Committed: Wed Mar 12 15:12:24 2003 UTC (21 years, 5 months ago) by mmeineke
File size: 12315 byte(s)
Log Message:
added the force parameter files into the distribution

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