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 1638 by chrisfen, Fri Oct 22 23:00:06 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/UseTheForce/LJFF.cpp (file contents), Revision 1741 by tim, Tue Nov 16 02:07:14 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 <  
41 < #ifdef IS_MPI
42 <  
43 <  MPI_Datatype mpiAtomStructType;
44 <  
45 < #endif
53 >    while(ffStream.getline(line, bufferSize)){
54 >        ++lineNo;
55  
56 <  class LinkedAtomType {
57 <  public:
58 <    LinkedAtomType(){
59 <      next = NULL;
60 <      name[0] = '\0';
61 <    }
62 <    ~LinkedAtomType(){ if( next != NULL ) delete next; }
56 >        //a line begins with '#' or '!' is comment which is skipped
57 >        if (line[0] != '#' ||line[0] != '!') {
58 >            StringTokenizer tokenizer(line);
59 >            
60 >            if (tokenizer.countToken() >= 4) {
61 >                atomTypeName = tokenizer.nextToken();
62 >                mass = tokenizer.nextTokenAsFloat();
63 >                epsilon = tokenizer.nextTokenAsFloat();
64 >                sigma = tokenizer.nextTokenAsFloat();
65  
66 <    LinkedAtomType* find(char* key){
67 <      if( !strcmp(name, key) ) return this;
68 <      if( next != NULL ) return next->find(key);
69 <      return NULL;
70 <    }
71 <    
66 >                atomType = new AtomType();
67 >                atomType->setName(atomTypeName);
68 >                atomType->setMass(mass);
69 >                
70 >                //by default, all of the properties in AtomTypeProperties is set to 0
71 >                //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
72 >                atomType->setLennardJones();
73  
74 <    void add( atomStruct &info ){
75 <    
76 <      // check for duplicates
77 <      
78 <      if( !strcmp( info.name, name ) ){
79 <        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 <    
86 <
87 < #ifdef IS_MPI
88 <    
89 <    void duplicate( atomStruct &info ){
90 <      strcpy(info.name, name);
91 <      info.mass     = mass;
92 <      info.epslon   = epslon;
93 <      info.sigma    = sigma;
94 <      info.ident    = ident;
95 <      info.last     = 0;
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" );
74 >                atomType->setIdent(ident);
75 >                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
76 >                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
77 >                atomType->complete();
78 >                //notify a new LJtype atom type is created
79 >                newLJtype(&ident, &sigma, &epsilon, &status);
80  
81 <    // we are now at the AtomTypes section.
82 <    
83 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
84 <    lineNum++;
85 <    
86 <    
87 <    // read a line, and start parseing out the atom types
81 >                //add atom type to AtomTypeContainer
82 >                addAtomType(atomTypeName, atomType);
83 >                ++ident;
84 >                
85 >            } else {
86 >                sprintf( painCave.errMsg,
87 >                       "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
88 >                       "in line %d : %s\n",
89 >                       filename.c_str(), lineNo, line);
90 >                painCave.severity = OOPSE_ERROR;
91 >                painCave.isFatal = 1;
92 >                simError();                
93 >            }
94 >            
95 >        }
96  
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();
97      }
98      
99 <    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 <      printf("Did setName\n");
362 <      at->setLennardJones();
363 <      printf("Did setLennardJones\n");
364 <      at->complete();
365 <      printf("Did complete\n");
366 <      
367 <    }
368 <    currentAtomType = currentAtomType->next;
369 <  }
370 <
371 <  currentAtomType = headAtomType;
372 <  while( currentAtomType != NULL ){
373 <    
374 <    if( currentAtomType->name[0] != '\0' ){
375 <      isError = 0;
376 <      newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
377 <                 &(currentAtomType->epslon), &isError);
378 <      if( isError ){
379 <        sprintf( painCave.errMsg,
380 <                 "Error initializing the \"%s\" LJ type in fortran\n",
381 <                 currentAtomType->name );
382 <        painCave.isFatal = 1;
383 <        simError();
384 <      }
385 <    }
386 <    
387 <    currentAtomType = currentAtomType->next;
388 <  }
389 <  
390 <  entry_plug->useLennardJones = 1;
391 <
392 < #ifdef IS_MPI
393 <  sprintf( checkPointMsg,
394 <           "LJFF atom structures successfully sent to fortran\n" );
395 <  MPIcheckPoint();
396 < #endif // is_mpi
397 <
99 >    delete ffStream;
100   }
101  
400 double LJFF::getAtomTypeMass (char* atomType) {
102  
103 <  currentAtomType = headAtomType->find( atomType );
403 <  if( currentAtomType == NULL ){
404 <    sprintf( painCave.errMsg,
405 <            "AtomType error, %s not found in force file.\n",
406 <             atomType );
407 <    painCave.isFatal = 1;
408 <    simError();
409 <  }
410 <
411 <  return currentAtomType->mass;
412 < }
413 <
414 < void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
415 <  
416 <  int i;
417 <
418 <  // initialize the atoms
419 <  
420 <
421 <  for( i=0; i<nAtoms; i++ ){
422 <    
423 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
424 <    if( currentAtomType == NULL ){
425 <      sprintf( painCave.errMsg,
426 <               "AtomType error, %s not found in force file.\n",
427 <               the_atoms[i]->getType() );
428 <      painCave.isFatal = 1;
429 <      simError();
430 <    }
431 <    
432 <    the_atoms[i]->setMass( currentAtomType->mass );
433 <    the_atoms[i]->setIdent( currentAtomType->ident );
434 <
435 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
436 <  }
437 < }
438 <
439 < void LJFF::initializeBonds( int nBonds, Bond** BondArray,
440 <                             bond_pair* the_bonds ){
441 <  
442 <    if( nBonds ){
443 <      sprintf( painCave.errMsg,
444 <               "LJFF does not support bonds.\n" );
445 <      painCave.isFatal = 1;
446 <      simError();
447 <    }
448 < }
449 <
450 < void LJFF::initializeBends( int nBends, Bend** bendArray,
451 <                             bend_set* the_bends ){
452 <
453 <    if( nBends ){
454 <      sprintf( painCave.errMsg,
455 <               "LJFF does not support bends.\n" );
456 <      painCave.isFatal = 1;
457 <      simError();
458 <    }
459 < }
460 <
461 < void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
462 <                                torsion_set* the_torsions ){
463 <
464 <    if( nTorsions ){
465 <      sprintf( painCave.errMsg,
466 <               "LJFF does not support torsions.\n" );
467 <      painCave.isFatal = 1;
468 <      simError();
469 <    }
470 < }
471 <
472 < void LJFF::fastForward( char* stopText, char* searchOwner ){
473 <
474 <  int foundText = 0;
475 <  char* the_token;
476 <
477 <  rewind( frcFile );
478 <  lineNum = 0;
479 <
480 <  eof_test = fgets( readLine, sizeof(readLine), frcFile );
481 <  lineNum++;
482 <  if( eof_test == NULL ){
483 <    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
484 <             " file is empty.\n",
485 <             searchOwner );
486 <    painCave.isFatal = 1;
487 <    simError();
488 <  }
489 <  
490 <  
491 <  while( !foundText ){
492 <    while( eof_test != NULL && readLine[0] != '#' ){
493 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
494 <      lineNum++;
495 <    }
496 <    if( eof_test == NULL ){
497 <      sprintf( painCave.errMsg,
498 <               "Error fast forwarding force file for %s at "
499 <               "line %d: file ended unexpectedly.\n",
500 <               searchOwner,
501 <               lineNum );
502 <      painCave.isFatal = 1;
503 <      simError();
504 <    }
505 <    
506 <    the_token = strtok( readLine, " ,;\t#\n" );
507 <    foundText = !strcmp( stopText, the_token );
508 <    
509 <    if( !foundText ){
510 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
511 <      lineNum++;
512 <      
513 <      if( eof_test == NULL ){
514 <        sprintf( painCave.errMsg,
515 <                 "Error fast forwarding force file for %s at "
516 <                 "line %d: file ended unexpectedly.\n",
517 <                 searchOwner,
518 <                 lineNum );
519 <        painCave.isFatal = 1;
520 <        simError();
521 <      }
522 <    }
523 <  }  
524 < }
525 <
526 <
527 <
528 < int LJ_NS::parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){
529 <
530 <  char* the_token;
531 <  
532 <  the_token = strtok( lineBuffer, " \n\t,;" );
533 <  if( the_token != NULL ){
534 <    
535 <    strcpy( info.name, the_token );
536 <          
537 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
538 <      sprintf( painCave.errMsg,
539 <               "Error parseing AtomTypes: line %d\n", lineNum );
540 <      painCave.isFatal = 1;
541 <      simError();
542 <    }
543 <    
544 <    info.mass = atof( the_token );
545 <    
546 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
547 <      sprintf( painCave.errMsg,
548 <               "Error parseing AtomTypes: line %d\n", lineNum );
549 <      painCave.isFatal = 1;
550 <      simError();
551 <    }
552 <        
553 <    info.epslon = atof( the_token );
554 <          
555 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
556 <      sprintf( painCave.errMsg,
557 <               "Error parseing AtomTypes: line %d\n", lineNum );
558 <      painCave.isFatal = 1;
559 <      simError();
560 <    }
561 <        
562 <    info.sigma = atof( the_token );
563 <    
564 <    return 1;
565 <  }
566 <  else return 0;
567 < }
568 <
569 <
103 > } //end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines