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:
trunk/OOPSE-3.0/src/UseTheForce/LJFF.cpp (file contents), Revision 1636 by chrisfen, Fri Oct 22 22:54:01 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/UseTheForce/LJFF.cpp (file contents), 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(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 <    }
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 <
99 < #endif
100 <
101 <    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(){
122 <
123 <  char fileName[200];
124 <  char* ffPath_env = "FORCE_PARAM_PATH";
125 <  char* ffPath;
126 <  char temp[200];
127 <
128 <  headAtomType = NULL;
129 <  currentAtomType = NULL;
130 <
131 <  // do the funtion wrapping
132 < // wrapMeFF( this );
133 <
134 < #ifdef IS_MPI
135 <  int i;
136 <  
137 <   // **********************************************************************
138 <  // Init the atomStruct mpi type
139 <
140 <  atomStruct atomProto; // mpiPrototype
141 <  int atomBC[3] = {15,3,2};  // block counts
142 <  MPI_Aint atomDspls[3];           // displacements
143 <  MPI_Datatype atomMbrTypes[3];    // member mpi types
144 <
145 <  MPI_Address(&atomProto.name, &atomDspls[0]);
146 <  MPI_Address(&atomProto.mass, &atomDspls[1]);
147 <  MPI_Address(&atomProto.ident, &atomDspls[2]);
148 <  
149 <  atomMbrTypes[0] = MPI_CHAR;
150 <  atomMbrTypes[1] = MPI_DOUBLE;
151 <  atomMbrTypes[2] = MPI_INT;
152 <  
153 <  for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
154 <  
155 <  MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
156 <  MPI_Type_commit(&mpiAtomStructType);
157 <
158 <  // ***********************************************************************
159 <  
160 <  if( worldRank == 0 ){
161 < #endif
162 <    
163 <    // generate the force file name
164 <    
165 <    strcpy( fileName, "LJFF.frc" );
166 <    // fprintf( stderr,"Trying to open %s\n", fileName );
167 <    
168 <    // attempt to open the file in the current directory first.
169 <    
170 <    frcFile = fopen( fileName, "r" );
171 <    
172 <    if( frcFile == NULL ){
173 <      
174 <      // next see if the force path enviorment variable is set
175 <      
176 <      ffPath = getenv( ffPath_env );
177 <      if( ffPath == NULL ) {
178 <        STR_DEFINE(ffPath, FRC_PATH );
179 <      }
180 <      
181 <      
182 <      strcpy( temp, ffPath );
183 <      strcat( temp, "/" );
184 <      strcat( temp, fileName );
185 <      strcpy( fileName, temp );
186 <      
187 <      frcFile = fopen( fileName, "r" );
188 <      
189 <      if( frcFile == NULL ){
190 <
191 <        sprintf( painCave.errMsg,
192 <                 "Error opening the force field parameter file:\n"
193 <                 "\t%s\n"
194 <                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
195 <                 "variable?\n",
196 <                 fileName );
197 <        painCave.severity = OOPSE_ERROR;
198 <        painCave.isFatal = 1;
199 <        simError();
200 <      }
201 <    }
202 <    
203 < #ifdef IS_MPI
204 <  }
205 <  
206 <  sprintf( checkPointMsg, "LJFF file opened sucessfully." );
207 <  MPIcheckPoint();
208 <  
209 < #endif // is_mpi
210 < }
211 <
212 <
213 < LJFF::~LJFF(){
214 <
215 <  if( headAtomType != NULL ) delete headAtomType;
216 <
217 < #ifdef IS_MPI
218 <  if( worldRank == 0 ){
219 < #endif // is_mpi
220 <    
221 <    fclose( frcFile );
222 <    
223 < #ifdef IS_MPI
224 <  }
225 < #endif // is_mpi
226 < }
227 <
228 < void LJFF::initForceField(){
229 <  
230 <  initFortran(0);
231 < }
232 <
233 < void LJFF::cleanMe( void ){
234 <
235 < #ifdef IS_MPI
236 <  
237 <  // keep the linked list in the mpi version
238 <
239 < #else // is_mpi
240 <
241 <  // delete the linked list in the single processor version
242 <
243 <  if( headAtomType != NULL ) delete headAtomType;
244 <
245 < #endif // is_mpi
246 < }
247 <
248 < void LJFF::readParams( void ){
249 <
250 <  atomStruct info;
251 <  info.last = 1; // initialize last to have the last set.
252 <                 // if things go well, last will be set to 0
253 <
254 <  int identNum;
255 <  int isError;
256 <
257 <  bigSigma = 0.0;
258 < #ifdef IS_MPI
259 <  if( worldRank == 0 ){
260 < #endif
261 <    
262 <    // read in the atom types.
263 <
264 <    headAtomType = new LinkedAtomType;
265 <    
266 <    fastForward( "AtomTypes", "initializeAtoms" );
267 <
268 <    // we are now at the AtomTypes section.
269 <    
270 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
271 <    lineNum++;
272 <    
273 <    
274 <    // read a line, and start parseing out the atom types
275 <
276 <    if( eof_test == NULL ){
277 <      sprintf( painCave.errMsg,
278 <               "Error in reading Atoms from force file at line %d.\n",
279 <               lineNum );
280 <      painCave.isFatal = 1;
281 <      simError();
98 >         ffStream.getline(line, bufferSize);
99      }
100      
101 <    identNum = 1;
285 <    // stop reading at end of file, or at next section
286 <    while( readLine[0] != '#' && eof_test != NULL ){
287 <
288 <      // toss comment lines
289 <      if( readLine[0] != '!' ){
290 <        
291 <        // the parser returns 0 if the line was blank
292 <        if( parseAtom( readLine, lineNum, info ) ){
293 <          info.ident = identNum;
294 <          headAtomType->add( info );;
295 <          identNum++;
296 <        }
297 <      }
298 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
299 <      lineNum++;
300 <    }
301 <
302 < #ifdef IS_MPI
303 <    
304 <    // send out the linked list to all the other processes
305 <
306 <    sprintf( checkPointMsg,
307 <             "LJFF atom structures read successfully." );
308 <    MPIcheckPoint();
309 <
310 <    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
311 <    while( currentAtomType != NULL ){
312 <      currentAtomType->duplicate( info );
313 <
314 <
315 <
316 <      sendFrcStruct( &info, mpiAtomStructType );
317 <
318 <      sprintf( checkPointMsg,
319 <               "successfully sent lJ force type: \"%s\"\n",
320 <               info.name );
321 <      MPIcheckPoint();
322 <
323 <      currentAtomType = currentAtomType->next;
324 <    }
325 <    info.last = 1;
326 <    sendFrcStruct( &info, mpiAtomStructType );
327 <    
328 <  }
329 <
330 <  else{
331 <    
332 <    // listen for node 0 to send out the force params
333 <    
334 <    MPIcheckPoint();
335 <
336 <    headAtomType = new LinkedAtomType;
337 <    receiveFrcStruct( &info, mpiAtomStructType );
338 <    
339 <    while( !info.last ){
340 <
341 <
342 <
343 <      headAtomType->add( info );
344 <      
345 <      MPIcheckPoint();
346 <
347 <      receiveFrcStruct( &info, mpiAtomStructType );
348 <    }
349 <  }
350 < #endif // is_mpi
351 <  
352 <  currentAtomType = headAtomType;
353 <  while( currentAtomType != NULL ){
354 <    
355 <    if( currentAtomType->name[0] != '\0' ){
356 <
357 <      AtomType* at = new AtomType();
358 <      at->setIdent(currentAtomType->ident);
359 <      printf ("currentName = %s\n", currentAtomType->name);
360 <      at->setName(currentAtomType->name);    
361 <      at->setLennardJones();
362 <      at->complete();
363 <      
364 <    }
365 <    currentAtomType = currentAtomType->next;
366 <  }
367 <
368 <  currentAtomType = headAtomType;
369 <  while( currentAtomType != NULL ){
370 <    
371 <    if( currentAtomType->name[0] != '\0' ){
372 <      isError = 0;
373 <      newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
374 <                 &(currentAtomType->epslon), &isError);
375 <      if( isError ){
376 <        sprintf( painCave.errMsg,
377 <                 "Error initializing the \"%s\" LJ type in fortran\n",
378 <                 currentAtomType->name );
379 <        painCave.isFatal = 1;
380 <        simError();
381 <      }
382 <    }
383 <    
384 <    currentAtomType = currentAtomType->next;
385 <  }
386 <  
387 <  entry_plug->useLennardJones = 1;
388 <
389 < #ifdef IS_MPI
390 <  sprintf( checkPointMsg,
391 <           "LJFF atom structures successfully sent to fortran\n" );
392 <  MPIcheckPoint();
393 < #endif // is_mpi
394 <
101 >    delete ffStream;
102   }
103  
397 double LJFF::getAtomTypeMass (char* atomType) {
104  
105 <  currentAtomType = headAtomType->find( atomType );
400 <  if( currentAtomType == NULL ){
401 <    sprintf( painCave.errMsg,
402 <            "AtomType error, %s not found in force file.\n",
403 <             atomType );
404 <    painCave.isFatal = 1;
405 <    simError();
406 <  }
407 <
408 <  return currentAtomType->mass;
409 < }
410 <
411 < void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
412 <  
413 <  int i;
414 <
415 <  // initialize the atoms
416 <  
417 <
418 <  for( i=0; i<nAtoms; i++ ){
419 <    
420 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
421 <    if( currentAtomType == NULL ){
422 <      sprintf( painCave.errMsg,
423 <               "AtomType error, %s not found in force file.\n",
424 <               the_atoms[i]->getType() );
425 <      painCave.isFatal = 1;
426 <      simError();
427 <    }
428 <    
429 <    the_atoms[i]->setMass( currentAtomType->mass );
430 <    the_atoms[i]->setIdent( currentAtomType->ident );
431 <
432 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
433 <  }
434 < }
435 <
436 < void LJFF::initializeBonds( int nBonds, Bond** BondArray,
437 <                             bond_pair* the_bonds ){
438 <  
439 <    if( nBonds ){
440 <      sprintf( painCave.errMsg,
441 <               "LJFF does not support bonds.\n" );
442 <      painCave.isFatal = 1;
443 <      simError();
444 <    }
445 < }
446 <
447 < void LJFF::initializeBends( int nBends, Bend** bendArray,
448 <                             bend_set* the_bends ){
449 <
450 <    if( nBends ){
451 <      sprintf( painCave.errMsg,
452 <               "LJFF does not support bends.\n" );
453 <      painCave.isFatal = 1;
454 <      simError();
455 <    }
456 < }
457 <
458 < void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
459 <                                torsion_set* the_torsions ){
460 <
461 <    if( nTorsions ){
462 <      sprintf( painCave.errMsg,
463 <               "LJFF does not support torsions.\n" );
464 <      painCave.isFatal = 1;
465 <      simError();
466 <    }
467 < }
468 <
469 < void LJFF::fastForward( char* stopText, char* searchOwner ){
470 <
471 <  int foundText = 0;
472 <  char* the_token;
473 <
474 <  rewind( frcFile );
475 <  lineNum = 0;
476 <
477 <  eof_test = fgets( readLine, sizeof(readLine), frcFile );
478 <  lineNum++;
479 <  if( eof_test == NULL ){
480 <    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
481 <             " file is empty.\n",
482 <             searchOwner );
483 <    painCave.isFatal = 1;
484 <    simError();
485 <  }
486 <  
487 <  
488 <  while( !foundText ){
489 <    while( eof_test != NULL && readLine[0] != '#' ){
490 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
491 <      lineNum++;
492 <    }
493 <    if( eof_test == NULL ){
494 <      sprintf( painCave.errMsg,
495 <               "Error fast forwarding force file for %s at "
496 <               "line %d: file ended unexpectedly.\n",
497 <               searchOwner,
498 <               lineNum );
499 <      painCave.isFatal = 1;
500 <      simError();
501 <    }
502 <    
503 <    the_token = strtok( readLine, " ,;\t#\n" );
504 <    foundText = !strcmp( stopText, the_token );
505 <    
506 <    if( !foundText ){
507 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
508 <      lineNum++;
509 <      
510 <      if( eof_test == NULL ){
511 <        sprintf( painCave.errMsg,
512 <                 "Error fast forwarding force file for %s at "
513 <                 "line %d: file ended unexpectedly.\n",
514 <                 searchOwner,
515 <                 lineNum );
516 <        painCave.isFatal = 1;
517 <        simError();
518 <      }
519 <    }
520 <  }  
521 < }
522 <
523 <
524 <
525 < int LJ_NS::parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){
526 <
527 <  char* the_token;
528 <  
529 <  the_token = strtok( lineBuffer, " \n\t,;" );
530 <  if( the_token != NULL ){
531 <    
532 <    strcpy( info.name, the_token );
533 <          
534 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
535 <      sprintf( painCave.errMsg,
536 <               "Error parseing AtomTypes: line %d\n", lineNum );
537 <      painCave.isFatal = 1;
538 <      simError();
539 <    }
540 <    
541 <    info.mass = atof( the_token );
542 <    
543 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
544 <      sprintf( painCave.errMsg,
545 <               "Error parseing AtomTypes: line %d\n", lineNum );
546 <      painCave.isFatal = 1;
547 <      simError();
548 <    }
549 <        
550 <    info.epslon = atof( the_token );
551 <          
552 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
553 <      sprintf( painCave.errMsg,
554 <               "Error parseing AtomTypes: line %d\n", lineNum );
555 <      painCave.isFatal = 1;
556 <      simError();
557 <    }
558 <        
559 <    info.sigma = atof( the_token );
560 <    
561 <    return 1;
562 <  }
563 <  else return 0;
564 < }
565 <
566 <
105 > } //end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines