ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/UseTheForce/LJFF.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/UseTheForce/LJFF.cpp (file contents):
Revision 1739 by tim, Wed Nov 3 18:00:36 2004 UTC vs.
Revision 1740 by tim, Mon Nov 15 23:00:32 2004 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
2 < #include <stdio.h>
3 < #include <string.h>
1 > /*
2 > * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 > *
4 > * Contact: oopse@oopse.org
5 > *
6 > * This program is free software; you can redistribute it and/or
7 > * modify it under the terms of the GNU Lesser General Public License
8 > * as published by the Free Software Foundation; either version 2.1
9 > * of the License, or (at your option) any later version.
10 > * All we ask is that proper credit is given for our work, which includes
11 > * - but is not limited to - adding the above copyright notice to the beginning
12 > * of your source code files, and to any copyright notice that you may distribute
13 > * with programs based on this work.
14 > *
15 > * This program is distributed in the hope that it will be useful,
16 > * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 > * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 > * GNU Lesser General Public License for more details.
19 > *
20 > * You should have received a copy of the GNU Lesser General Public License
21 > * along with this program; if not, write to the Free Software
22 > * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23 > *
24 > */
25  
26 < #include <iostream>
6 < using namespace std;
7 <
8 < #ifdef IS_MPI
9 < #include <mpi.h>
10 < #endif //is_mpi
11 <
12 < #include "UseTheForce/ForceFields.hpp"
13 < #include "primitives/SRI.hpp"
14 < #include "utils/simError.h"
15 < #include "types/AtomType.hpp"
26 > #include "UseTheForce/LJFF.hpp"
27   #include "UseTheForce/DarkSide/lj_interface.h"
28  
29 + namespace oopse {
30  
31 < #ifdef IS_MPI
32 < #include "UseTheForce/mpiForceField.h"
33 < #endif // is_mpi
31 > //definition of createLJFF
32 > ForceField* createLJFF() {
33 >    return new LJFF();
34 > }
35  
36 + //register createLJFF to ForceFieldCreator
37 + ForceFieldFactory::getInstance()->registerForceField("LJFF", createLJFF);
38  
39 <
40 < namespace LJ_NS{
41 <
42 <  // Declare the structures that will be passed by the parser and  MPI
43 <  
44 <  typedef struct{
45 <    char name[15];
39 > void LJFF::parse(const std::string& filename) {
40 >    ifstrstream* ffStream;
41 >    ffStream = openForceFiledFile(filename);
42 >    const int bufferSize = 65535;
43 >    char line[bufferSize];
44 >    AtomType* atomType;
45 >    std::string atomTypeName;
46      double mass;
47 <    double epslon;
47 >    double epsilon;
48      double sigma;
49 <    int ident;
50 <    int last;      //  0  -> default
51 <                   //  1  -> in MPI: tells nodes to stop listening
37 <  } atomStruct;
49 >    int status;
50 >    int ident = 1; //ident begins from 1 (c++ does not need to know it, fortran's index begins from 1)
51 >    int lineNo = 0;
52  
53 <  int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
54 <  
55 < #ifdef IS_MPI
42 <  
43 <  MPI_Datatype mpiAtomStructType;
44 <  
45 < #endif
53 >    ffStream.getline(line, bufferSize);    
54 >    while( !ffStream.eof() ){
55 >        ++lineNo;
56  
57 <  class LinkedAtomType {
58 <  public:
59 <    LinkedAtomType(){
60 <      next = NULL;
61 <      name[0] = '\0';
62 <    }
63 <    ~LinkedAtomType(){ if( next != NULL ) delete next; }
57 >        //a line begins with '#' or '!' is comment which is skipped
58 >        if (line[0] != '#' ||line[0] != '!') {
59 >            StringTokenizer tokenizer(line);
60 >            
61 >            if (tokenizer.countToken() >= 4) {
62 >                atomTypeName = tokenizer.nextToken();
63 >                mass = tokenizer.nextTokenAsFloat();
64 >                epsilon = tokenizer.nextTokenAsFloat();
65 >                sigma = tokenizer.nextTokenAsFloat();
66  
67 <    LinkedAtomType* find(const char* key){
68 <      if( !strcmp(name, key) ) return this;
69 <      if( next != NULL ) return next->find(key);
70 <      return NULL;
71 <    }
72 <    
67 >                atomType = new AtomType();
68 >                atomType->setName(atomTypeName);
69 >                atomType->setMass(mass);
70 >                
71 >                //by default, all of the properties in AtomTypeProperties is set to 0
72 >                //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
73 >                atomType->setLennardJones();
74  
75 <    void add( atomStruct &info ){
76 <    
77 <      // check for duplicates
78 <      
79 <      if( !strcmp( info.name, name ) ){
80 <        sprintf( painCave.errMsg,
68 <                 "Duplicate LJ atom type \"%s\" found in "
69 <                 "the LJFF param file./n",
70 <                 name );
71 <        painCave.isFatal = 1;
72 <        simError();
73 <      }
74 <      
75 <      if( next != NULL ) next->add(info);
76 <      else{
77 <        next = new LinkedAtomType();
78 <        strcpy(next->name, info.name);
79 <        next->mass     = info.mass;
80 <        next->epslon   = info.epslon;
81 <        next->sigma    = info.sigma;
82 <        next->ident    = info.ident;
83 <      }
84 <    }
85 <    
75 >                atomType->setIdent(ident);
76 >                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
77 >                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
78 >                atomType->complete();
79 >                //notify a new LJtype atom type is created
80 >                newLJtype(&ident, &sigma, &epsilon, &status);
81  
82 < #ifdef IS_MPI
83 <    
84 <    void duplicate( atomStruct &info ){
85 <      strcpy(info.name, name);
86 <      info.mass     = mass;
87 <      info.epslon   = epslon;
88 <      info.sigma    = sigma;
89 <      info.ident    = ident;
90 <      info.last     = 0;
91 <    }
92 <
93 <
94 < #endif
95 <
96 <    char name[15];
102 <    double mass;
103 <    double epslon;
104 <    double sigma;
105 <    int ident;
106 <    LinkedAtomType* next;
107 <  };
108 <
109 <  LinkedAtomType* headAtomType;
110 <  LinkedAtomType* currentAtomType;
111 <
112 < }
113 <
114 < using namespace LJ_NS;
115 <
116 < //****************************************************************
117 < // begins the actual forcefield stuff.  
118 < //****************************************************************
119 <
120 <
121 < LJFF::LJFF() : ForceFields() {
122 <
123 <  string fileName;
124 <  string tempString;
125 <    
126 <  headAtomType = NULL;
127 <  currentAtomType = NULL;
128 <
129 <  // do the funtion wrapping
130 < // wrapMeFF( this );
131 <
132 < #ifdef IS_MPI
133 <  int i;
134 <  
135 <   // **********************************************************************
136 <  // Init the atomStruct mpi type
137 <
138 <  atomStruct atomProto; // mpiPrototype
139 <  int atomBC[3] = {15,3,2};  // block counts
140 <  MPI_Aint atomDspls[3];           // displacements
141 <  MPI_Datatype atomMbrTypes[3];    // member mpi types
142 <
143 <  MPI_Address(&atomProto.name, &atomDspls[0]);
144 <  MPI_Address(&atomProto.mass, &atomDspls[1]);
145 <  MPI_Address(&atomProto.ident, &atomDspls[2]);
146 <  
147 <  atomMbrTypes[0] = MPI_CHAR;
148 <  atomMbrTypes[1] = MPI_DOUBLE;
149 <  atomMbrTypes[2] = MPI_INT;
150 <  
151 <  for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
152 <  
153 <  MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
154 <  MPI_Type_commit(&mpiAtomStructType);
155 <
156 <  // ***********************************************************************
157 <  
158 <  if( worldRank == 0 ){
159 < #endif
160 <    
161 <    // generate the force file name  
162 <
163 <    fileName = "LJFF.frc";
164 <    // fprintf( stderr,"Trying to open %s\n", fileName );
165 <    
166 <    // attempt to open the file in the current directory first.
167 <    
168 <    frcFile = fopen( fileName.c_str(), "r" );
169 <    
170 <    if( frcFile == NULL ){
171 <      
172 <      // next see if the force path enviorment variable is set
173 <
174 <      tempString = ffPath + "/" + fileName;
175 <      fileName = tempString;
176 <      
177 <      frcFile = fopen( fileName.c_str(), "r" );
178 <      
179 <      if( frcFile == NULL ){
180 <
181 <        sprintf( painCave.errMsg,
182 <                 "Error opening the force field parameter file:\n"
183 <                 "\t%s\n"
184 <                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
185 <                 "variable?\n",
186 <                 fileName.c_str() );
187 <        painCave.severity = OOPSE_ERROR;
188 <        painCave.isFatal = 1;
189 <        simError();
190 <      }
191 <    }
192 <    
193 < #ifdef IS_MPI
194 <  }
195 <  
196 <  sprintf( checkPointMsg, "LJFF file opened sucessfully." );
197 <  MPIcheckPoint();
198 <  
199 < #endif // is_mpi
200 < }
201 <
202 <
203 < LJFF::~LJFF(){
204 <
205 <  if( headAtomType != NULL ) delete headAtomType;
206 <
207 < #ifdef IS_MPI
208 <  if( worldRank == 0 ){
209 < #endif // is_mpi
210 <    
211 <    fclose( frcFile );
212 <    
213 < #ifdef IS_MPI
214 <  }
215 < #endif // is_mpi
216 < }
217 <
218 < void LJFF::initForceField(){
219 <  
220 <  initFortran(0);
221 < }
222 <
223 < void LJFF::cleanMe( void ){
224 <
225 < #ifdef IS_MPI
226 <  
227 <  // keep the linked list in the mpi version
228 <
229 < #else // is_mpi
230 <
231 <  // delete the linked list in the single processor version
232 <
233 <  if( headAtomType != NULL ) delete headAtomType;
234 <
235 < #endif // is_mpi
236 < }
237 <
238 < void LJFF::readParams( void ){
82 >                //add atom type to AtomTypeContainer
83 >                addAtomType(atomTypeName, atomType);
84 >                ++ident;
85 >                
86 >            } else {
87 >                sprintf( painCave.errMsg,
88 >                       "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
89 >                       "in line %d : %s\n",
90 >                       filename.c_str(), lineNo, line);
91 >                painCave.severity = OOPSE_ERROR;
92 >                painCave.isFatal = 1;
93 >                simError();                
94 >            }
95 >            
96 >        }
97  
98 <  atomStruct info;
241 <  info.last = 1; // initialize last to have the last set.
242 <                 // if things go well, last will be set to 0
243 <
244 <  int identNum;
245 <  int isError;
246 <
247 <  bigSigma = 0.0;
248 < #ifdef IS_MPI
249 <  if( worldRank == 0 ){
250 < #endif
251 <    
252 <    // read in the atom types.
253 <
254 <    headAtomType = new LinkedAtomType;
255 <    
256 <    fastForward( "AtomTypes", "initializeAtoms" );
257 <
258 <    // we are now at the AtomTypes section.
259 <    
260 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
261 <    lineNum++;
262 <    
263 <    
264 <    // read a line, and start parseing out the atom types
265 <
266 <    if( eof_test == NULL ){
267 <      sprintf( painCave.errMsg,
268 <               "Error in reading Atoms from force file at line %d.\n",
269 <               lineNum );
270 <      painCave.isFatal = 1;
271 <      simError();
98 >         ffStream.getline(line, bufferSize);
99      }
100      
101 <    identNum = 1;
275 <    // stop reading at end of file, or at next section
276 <    while( readLine[0] != '#' && eof_test != NULL ){
277 <
278 <      // toss comment lines
279 <      if( readLine[0] != '!' ){
280 <        
281 <        // the parser returns 0 if the line was blank
282 <        if( parseAtom( readLine, lineNum, info ) ){
283 <          info.ident = identNum;
284 <          headAtomType->add( info );;
285 <          identNum++;
286 <        }
287 <      }
288 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
289 <      lineNum++;
290 <    }
291 <
292 < #ifdef IS_MPI
293 <    
294 <    // send out the linked list to all the other processes
295 <
296 <    sprintf( checkPointMsg,
297 <             "LJFF atom structures read successfully." );
298 <    MPIcheckPoint();
299 <
300 <    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
301 <    while( currentAtomType != NULL ){
302 <      currentAtomType->duplicate( info );
303 <
304 <
305 <
306 <      sendFrcStruct( &info, mpiAtomStructType );
307 <
308 <      sprintf( checkPointMsg,
309 <               "successfully sent lJ force type: \"%s\"\n",
310 <               info.name );
311 <      MPIcheckPoint();
312 <
313 <      currentAtomType = currentAtomType->next;
314 <    }
315 <    info.last = 1;
316 <    sendFrcStruct( &info, mpiAtomStructType );
317 <    
318 <  }
319 <
320 <  else{
321 <    
322 <    // listen for node 0 to send out the force params
323 <    
324 <    MPIcheckPoint();
325 <
326 <    headAtomType = new LinkedAtomType;
327 <    receiveFrcStruct( &info, mpiAtomStructType );
328 <    
329 <    while( !info.last ){
330 <
331 <
332 <
333 <      headAtomType->add( info );
334 <      
335 <      MPIcheckPoint();
336 <
337 <      receiveFrcStruct( &info, mpiAtomStructType );
338 <    }
339 <  }
340 < #endif // is_mpi
341 <  
342 <  currentAtomType = headAtomType;
343 <  while( currentAtomType != NULL ){
344 <    
345 <    if( currentAtomType->name[0] != '\0' ){
346 <
347 <      AtomType* at = new AtomType();
348 <      at->setIdent(currentAtomType->ident);
349 <      printf ("currentName = %s\n", currentAtomType->name);
350 <      at->setName(currentAtomType->name);    
351 <      printf("Did setName\n");
352 <      at->setLennardJones();
353 <      printf("Did setLennardJones\n");
354 <      at->complete();
355 <      printf("Did complete\n");
356 <      
357 <    }
358 <    currentAtomType = currentAtomType->next;
359 <  }
360 <
361 <  currentAtomType = headAtomType;
362 <  while( currentAtomType != NULL ){
363 <    
364 <    if( currentAtomType->name[0] != '\0' ){
365 <      isError = 0;
366 <      newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
367 <                 &(currentAtomType->epslon), &isError);
368 <      if( isError ){
369 <        sprintf( painCave.errMsg,
370 <                 "Error initializing the \"%s\" LJ type in fortran\n",
371 <                 currentAtomType->name );
372 <        painCave.isFatal = 1;
373 <        simError();
374 <      }
375 <    }
376 <    
377 <    currentAtomType = currentAtomType->next;
378 <  }
379 <  
380 <  entry_plug->useLennardJones = 1;
381 <
382 < #ifdef IS_MPI
383 <  sprintf( checkPointMsg,
384 <           "LJFF atom structures successfully sent to fortran\n" );
385 <  MPIcheckPoint();
386 < #endif // is_mpi
387 <
101 >    delete ffStream;
102   }
103  
390 void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
391  
392  int i;
104  
105 <  // initialize the atoms
395 <  
396 <
397 <  for( i=0; i<nAtoms; i++ ){
398 <    
399 <    currentAtomType = headAtomType->find(the_atoms[i]->getType().c_str() );
400 <    if( currentAtomType == NULL ){
401 <      sprintf( painCave.errMsg,
402 <               "AtomType error, %s not found in force file.\n",
403 <               the_atoms[i]->getType().c_str() );
404 <      painCave.isFatal = 1;
405 <      simError();
406 <    }
407 <    
408 <    the_atoms[i]->setMass( currentAtomType->mass );
409 <    the_atoms[i]->setIdent( currentAtomType->ident );
410 <
411 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
412 <  }
413 < }
414 <
415 < void LJFF::initializeBonds( int nBonds, Bond** BondArray,
416 <                             bond_pair* the_bonds ){
417 <  
418 <    if( nBonds ){
419 <      sprintf( painCave.errMsg,
420 <               "LJFF does not support bonds.\n" );
421 <      painCave.isFatal = 1;
422 <      simError();
423 <    }
424 < }
425 <
426 < void LJFF::initializeBends( int nBends, Bend** bendArray,
427 <                             bend_set* the_bends ){
428 <
429 <    if( nBends ){
430 <      sprintf( painCave.errMsg,
431 <               "LJFF does not support bends.\n" );
432 <      painCave.isFatal = 1;
433 <      simError();
434 <    }
435 < }
436 <
437 < void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
438 <                                torsion_set* the_torsions ){
439 <
440 <    if( nTorsions ){
441 <      sprintf( painCave.errMsg,
442 <               "LJFF does not support torsions.\n" );
443 <      painCave.isFatal = 1;
444 <      simError();
445 <    }
446 < }
447 <
448 < void LJFF::fastForward( char* stopText, char* searchOwner ){
449 <
450 <  int foundText = 0;
451 <  char* the_token;
452 <
453 <  rewind( frcFile );
454 <  lineNum = 0;
455 <
456 <  eof_test = fgets( readLine, sizeof(readLine), frcFile );
457 <  lineNum++;
458 <  if( eof_test == NULL ){
459 <    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
460 <             " file is empty.\n",
461 <             searchOwner );
462 <    painCave.isFatal = 1;
463 <    simError();
464 <  }
465 <  
466 <  
467 <  while( !foundText ){
468 <    while( eof_test != NULL && readLine[0] != '#' ){
469 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
470 <      lineNum++;
471 <    }
472 <    if( eof_test == NULL ){
473 <      sprintf( painCave.errMsg,
474 <               "Error fast forwarding force file for %s at "
475 <               "line %d: file ended unexpectedly.\n",
476 <               searchOwner,
477 <               lineNum );
478 <      painCave.isFatal = 1;
479 <      simError();
480 <    }
481 <    
482 <    the_token = strtok( readLine, " ,;\t#\n" );
483 <    foundText = !strcmp( stopText, the_token );
484 <    
485 <    if( !foundText ){
486 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
487 <      lineNum++;
488 <      
489 <      if( eof_test == NULL ){
490 <        sprintf( painCave.errMsg,
491 <                 "Error fast forwarding force file for %s at "
492 <                 "line %d: file ended unexpectedly.\n",
493 <                 searchOwner,
494 <                 lineNum );
495 <        painCave.isFatal = 1;
496 <        simError();
497 <      }
498 <    }
499 <  }  
500 < }
501 <
502 <
503 <
504 < int LJ_NS::parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){
505 <
506 <  char* the_token;
507 <  
508 <  the_token = strtok( lineBuffer, " \n\t,;" );
509 <  if( the_token != NULL ){
510 <    
511 <    strcpy( info.name, the_token );
512 <          
513 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
514 <      sprintf( painCave.errMsg,
515 <               "Error parseing AtomTypes: line %d\n", lineNum );
516 <      painCave.isFatal = 1;
517 <      simError();
518 <    }
519 <    
520 <    info.mass = atof( the_token );
521 <    
522 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
523 <      sprintf( painCave.errMsg,
524 <               "Error parseing AtomTypes: line %d\n", lineNum );
525 <      painCave.isFatal = 1;
526 <      simError();
527 <    }
528 <        
529 <    info.epslon = atof( the_token );
530 <          
531 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
532 <      sprintf( painCave.errMsg,
533 <               "Error parseing AtomTypes: line %d\n", lineNum );
534 <      painCave.isFatal = 1;
535 <      simError();
536 <    }
537 <        
538 <    info.sigma = atof( the_token );
539 <    
540 <    return 1;
541 <  }
542 <  else return 0;
543 < }
544 <
545 <
105 > } //end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines