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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines