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

Comparing trunk/OOPSE-3.0/src/UseTheForce/Shapes_FF.cpp (file contents):
Revision 1650 by gezelter, Tue Oct 26 22:24:52 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 < /*
2 < *  Shapes_FF.cpp
3 < *  oopse
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4 < *  Created by Chris Fennell on 10/20/04.
5 < *  Copyright 2004 __MyCompanyName__. All rights reserved.
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8   *
9 + * 1. Acknowledgement of the program authors must be made in any
10 + *    publication of scientific results based in part on use of the
11 + *    program.  An acceptable form of acknowledgement is citation of
12 + *    the article in which the program was described (Matthew
13 + *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + *    Parallel Simulation Engine for Molecular Dynamics,"
16 + *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + *
18 + * 2. Redistributions of source code must retain the above copyright
19 + *    notice, this list of conditions and the following disclaimer.
20 + *
21 + * 3. Redistributions in binary form must reproduce the above copyright
22 + *    notice, this list of conditions and the following disclaimer in the
23 + *    documentation and/or other materials provided with the
24 + *    distribution.
25 + *
26 + * This software is provided "AS IS," without a warranty of any
27 + * kind. All express or implied conditions, representations and
28 + * warranties, including any implied warranty of merchantability,
29 + * fitness for a particular purpose or non-infringement, are hereby
30 + * excluded.  The University of Notre Dame and its licensors shall not
31 + * be liable for any damages suffered by licensee as a result of
32 + * using, modifying or distributing the software or its
33 + * derivatives. In no event will the University of Notre Dame or its
34 + * licensors be liable for any lost revenue, profit or data, or for
35 + * direct, indirect, special, consequential, incidental or punitive
36 + * damages, however caused and regardless of the theory of liability,
37 + * arising out of the use of or inability to use software, even if the
38 + * University of Notre Dame has been advised of the possibility of
39 + * such damages.
40   */
41 <
41 >
42   #include <stdlib.h>
43   #include <stdio.h>
44   #include <string.h>
# Line 33 | Line 65 | using namespace std;
65   #include "UseTheForce/mpiForceField.h"
66   #endif // is_mpi
67  
68 < using namespace std;
68 >
69   using namespace oopse;
70  
71   Shapes_FF::~Shapes_FF(){
# Line 53 | Line 85 | void Shapes_FF::calcRcut( void ){
85   void Shapes_FF::calcRcut( void ){
86    
87   #ifdef IS_MPI
88 <  double tempShapesRcut = shapesRcut;
88 >  double tempShapesRcut = bigContact;
89    MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
90                   MPI_COMM_WORLD);
91 + #else
92 +  shapesRcut = bigContact;
93   #endif  //is_mpi
94    entry_plug->setDefaultRcut(shapesRcut);
95   }
# Line 70 | Line 104 | void Shapes_FF::readParams( void ){
104  
105    char readLine[1024];
106  
107 <  string fileName;
108 <  string shapeFileName;
109 <  string tempString;
107 >   std::string fileName;
108 >   std::string shapeFileName;
109 >   std::string tempString;
110  
111    char *nameToken;
112    char *delim = " ,;\t\n";
# Line 82 | Line 116 | void Shapes_FF::readParams( void ){
116    int nStrength = 0;
117    int myATID;
118    int isError;
119 <  string nameString;
119 >   std::string nameString;
120    AtomType* at;
121    DirectionalAtomType* dat;
122    ShapeAtomType* st;
123  
124 <  map<string, AtomType*>::iterator iter;
124 >   std::map<string, AtomType*>::iterator iter;
125  
126    // vectors for shape transfer to fortran
127 <  vector<RealSphericalHarmonic*> tempSHVector;
128 <  vector<int> contactL;
129 <  vector<int> contactM;
130 <  vector<int> contactFunc;
131 <  vector<double> contactCoeff;
132 <  vector<int> rangeL;
133 <  vector<int> rangeM;
134 <  vector<int> rangeFunc;
135 <  vector<double> rangeCoeff;
136 <  vector<int> strengthL;
137 <  vector<int> strengthM;
138 <  vector<int> strengthFunc;
139 <  vector<double> strengthCoeff;
127 >   std::vector<RealSphericalHarmonic*> tempSHVector;
128 >   std::vector<int> contactL;
129 >   std::vector<int> contactM;
130 >   std::vector<int> contactFunc;
131 >   std::vector<double> contactCoeff;
132 >   std::vector<int> rangeL;
133 >   std::vector<int> rangeM;
134 >   std::vector<int> rangeFunc;
135 >   std::vector<double> rangeCoeff;
136 >   std::vector<int> strengthL;
137 >   std::vector<int> strengthM;
138 >   std::vector<int> strengthFunc;
139 >   std::vector<double> strengthCoeff;
140    
141    // generate the force file name  
142    fileName = "Shapes.frc";
# Line 170 | Line 204 | void Shapes_FF::readParams( void ){
204            st = new ShapeAtomType();
205  
206            st->setName(nameString);
207 <          myATID = atomTypeMap.size();
207 >          myATID = atomTypeMap.size() + 1;
208            st->setIdent(myATID);
209            parseShapeFile(shapeFileName, st);
210            st->complete();
211            atomTypeMap.insert(make_pair(nameString, st));
212            
213          } else {
214 <          // atomType map already contained this string (i.e. it was
214 >          // atomType map already contained this  std::string (i.e. it was
215            // declared in a previous block, and we just need to add
216            // the shape-specific information for this AtomType:
217  
# Line 257 | Line 291 | void Shapes_FF::readParams( void ){
291          
292          isError = 0;
293          myATID = at->getIdent();
294 +
295          makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
296                     &contactCoeff[0],
297                     &nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
# Line 265 | Line 300 | void Shapes_FF::readParams( void ){
300                     &strengthFunc[0], &strengthCoeff[0],
301                     &myATID,
302                     &isError);
303 +
304          if( isError ){
305            sprintf( painCave.errMsg,
306                     "Error initializing the \"%s\" shape in fortran\n",
# Line 274 | Line 310 | void Shapes_FF::readParams( void ){
310          }
311        }
312      }
313 +  }
314 +  
315 +  isError = 0;
316 +  completeShapeFF(&isError);
317 +  if( isError ){
318 +    sprintf( painCave.errMsg,
319 +             "Error completing Shape FF in fortran\n");
320 +    painCave.isFatal = 1;
321 +    simError();
322    }
323    
324   #ifdef IS_MPI
# Line 291 | Line 336 | void Shapes_FF::initializeAtoms( int nAtoms, Atom** th
336   void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
337    
338    int i,j,k;
339 <  map<string, AtomType*>::iterator iter;
339 >   std::map<string, AtomType*>::iterator iter;
340  
341    // initialize the atoms
342    DirectionalAtom* dAtom;
343    AtomType* at;
344    DirectionalAtomType* dat;
345    ShapeAtomType* sat;
346 <  double sigma;
346 >  double longCutoff;
347    double ji[3];
348    double inertialMat[3][3];
349    Mat3x3d momInt;
350 <  string myTypeString;
350 >   std::string myTypeString;
351  
352 +  bigContact = 0.0;
353 +
354    for( i=0; i<nAtoms; i++ ){
355      
356 <    myTypeString = the_atoms[i]->getType();
356 >    myTypeString = the_atoms[i]->getType().c_str();
357      iter = atomTypeMap.find(myTypeString);
358  
359      if (iter == atomTypeMap.end()) {
360        sprintf( painCave.errMsg,
361                 "AtomType error, %s not found in force file.\n",
362 <               the_atoms[i]->getType() );
362 >               the_atoms[i]->getType().c_str() );
363        painCave.isFatal = 1;
364        simError();
365      } else {
# Line 325 | Line 372 | void Shapes_FF::initializeAtoms( int nAtoms, Atom** th
372        if ( at->isShape() ) {
373          
374          sat = (ShapeAtomType*)at;
375 <        sigma = findLargestContactDistance(sat);
376 <        if (sigma > bigSigma) bigSigma = sigma;
377 <
375 >        longCutoff = findCutoffDistance(sat);
376 >        if (longCutoff > bigContact) bigContact = longCutoff;  
377 >        cout << bigContact << " is the cutoff value\n";
378 >        
379 >        entry_plug->useShapes = 1;
380        }
381  
382        the_atoms[i]->setHasCharge(at->isCharge());
# Line 336 | Line 385 | void Shapes_FF::initializeAtoms( int nAtoms, Atom** th
385  
386          dat = (DirectionalAtomType*)at;
387          dAtom = (DirectionalAtom *) the_atoms[i];
388 <
388 >        
389          momInt = dat->getI();
390  
391          // zero out the moments of inertia matrix
# Line 400 | Line 449 | void Shapes_FF::parseShapeFile(string shapeFileName, S
449    int nTokens;
450    Mat3x3d momInert;
451    RealSphericalHarmonic* rsh;
452 <  vector<RealSphericalHarmonic*> functionVector;
452 >   std::vector<RealSphericalHarmonic*> functionVector;
453    ifstrstream shapeFile;
454 <  string tempString;
454 >   std::string tempString;
455  
456    shapeFile.open( shapeFileName.c_str() );
457    
# Line 450 | Line 499 | void Shapes_FF::parseShapeFile(string shapeFileName, S
499            painCave.severity = OOPSE_ERROR;
500            painCave.isFatal = 1;
501            simError();
502 <
502 >        } else {
503            token = strtok(inLine, delim);
504            token = strtok(NULL, delim);
505            st->setMass(atof(token));
# Line 464 | Line 513 | void Shapes_FF::parseShapeFile(string shapeFileName, S
513          }
514        }
515      }
516 +    shapeFile.getline(inLine, MAXLEN);
517    }
518  
519    // now grab the contact functions
# Line 476 | Line 526 | void Shapes_FF::parseShapeFile(string shapeFileName, S
526      if( inLine[0] != '!' && inLine[0] != '#' ){
527        // end marks section completion
528        if (isEndLine(inLine)) break;
479 –      
529        nTokens = countTokens(inLine, delim);
530        if (nTokens != 0) {
531          if (nTokens < 4) {
# Line 486 | Line 535 | void Shapes_FF::parseShapeFile(string shapeFileName, S
535            painCave.severity = OOPSE_ERROR;
536            painCave.isFatal = 1;
537            simError();
538 <          
538 >        } else {
539            // read in a spherical harmonic function
540            token = strtok(inLine, delim);
541            rsh = new RealSphericalHarmonic();
# Line 505 | Line 554 | void Shapes_FF::parseShapeFile(string shapeFileName, S
554          }
555        }
556      }
557 +    shapeFile.getline(inLine, MAXLEN);
558    }
559  
560    // pass contact functions to ShapeType
561 +
562    st->setContactFuncs(functionVector);
563  
564    // now grab the range functions
# Line 530 | Line 581 | void Shapes_FF::parseShapeFile(string shapeFileName, S
581            painCave.severity = OOPSE_ERROR;
582            painCave.isFatal = 1;
583            simError();
584 +        } else {
585            
586            // read in a spherical harmonic function
587            token = strtok(inLine, delim);
# Line 550 | Line 602 | void Shapes_FF::parseShapeFile(string shapeFileName, S
602          }
603        }
604      }
605 +    shapeFile.getline(inLine, MAXLEN);
606    }
607  
608    // pass range functions to ShapeType
# Line 575 | Line 628 | void Shapes_FF::parseShapeFile(string shapeFileName, S
628            painCave.severity = OOPSE_ERROR;
629            painCave.isFatal = 1;
630            simError();
631 +        } else {
632            
633            // read in a spherical harmonic function
634            token = strtok(inLine, delim);
# Line 594 | Line 648 | void Shapes_FF::parseShapeFile(string shapeFileName, S
648          }
649        }
650      }
651 +    shapeFile.getline(inLine, MAXLEN);
652    }
653  
654    // pass strength functions to ShapeType
# Line 607 | Line 662 | double Shapes_FF::findLargestContactDistance(ShapeAtom
662    int i, j,  nSteps;
663    double theta, thetaStep, thetaMin, costheta;
664    double phi, phiStep;
665 <  double sigma, bigSigma;
665 >  double sigma, bs;
666 >  
667 >  nSteps = 16;
668  
612 –  nSteps = 16;
613 –
669    thetaStep = M_PI / nSteps;
670    thetaMin = thetaStep / 2.0;
671    phiStep = thetaStep * 2.0;
672 <  bigSigma = 0.0;
672 >  bs = 0.0;
673    
674    for (i = 0; i < nSteps; i++) {
675      
# Line 627 | Line 682 | double Shapes_FF::findLargestContactDistance(ShapeAtom
682  
683        sigma = st->getContactValueAt(costheta, phi);
684        
685 <      if (sigma > bigSigma) bigSigma = sigma;
685 >      if (sigma > bs) bs = sigma;
686      }
687    }
688  
689 <  return bigSigma;  
689 >  return bs;  
690   }
691 +
692 +
693 + double Shapes_FF::findCutoffDistance(ShapeAtomType* st) {
694 +  int i, j,  nSteps;
695 +  double theta, thetaStep, thetaMin, costheta;
696 +  double phi, phiStep;
697 +  double sigma, range;
698 +  double bigCut, tempCut;
699 +
700 +  nSteps = 16;
701 +
702 +  thetaStep = M_PI / nSteps;
703 +  thetaMin = thetaStep / 2.0;
704 +  phiStep = thetaStep * 2.0;
705 +  bigCut = 0.0;
706 +  
707 +  for (i = 0; i < nSteps; i++) {
708 +    
709 +    theta = thetaMin + i * thetaStep;
710 +    costheta = cos(theta);
711 +
712 +    for (j = 0; j < nSteps; j++) {
713 +
714 +      phi = j*phiStep;
715 +
716 +      sigma = st->getContactValueAt(costheta, phi);
717 +      range = st->getRangeValueAt(costheta, phi);
718 +
719 +       // cutoff for a shape is taken to be (2.5*rangeVal + contactVal)
720 +      tempCut = 1.5*range + sigma;
721 +
722 +      if (tempCut > bigCut) bigCut = tempCut;
723 +    }
724 +  }
725 +
726 +  return bigCut;  
727 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines