ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/SHAPES/forcer.cpp
(Generate patch)

Comparing trunk/SHAPES/forcer.cpp (file contents):
Revision 1244 by gezelter, Fri Jun 4 16:21:23 2004 UTC vs.
Revision 1292 by chrisfen, Wed Jun 23 23:19:43 2004 UTC

# Line 4 | Line 4
4   #include <vector>
5   #include "forcerCmd.h"
6   #include "PDBReader.hpp"
7 + #include "RigidBody.hpp"
8 + #include "GridBuilder.hpp"
9 + #include "SphereHarm.hpp"
10  
11   #define MK_STR(s) # s
12   #define STR_DEFINE(t, s) t = MK_STR(s)
# Line 47 | Line 50 | int main(int argc, char* argv[]){
50    vector<Atype*>   vdwAtypes;
51    vector<Atype*>::iterator i;
52    Atype* at;
53 +  RigidBody* rb;
54    vector<VDWAtom*> theAtoms;
55    vector<VDWAtom*>::iterator j;
56    VDWAtom* atom;
57 <
57 >  GridBuilder* gb;
58 >  vector<double> sigmaGrid;
59 >  vector<double> epsGrid;
60 >  vector<double> sGrid;
61 >  SphereHarm* harmonize;
62 >  
63    double mass, rpar, eps;
64 +  double xyz3[3];
65 +  double moments[3][3];
66    string fileName;  
67    char vdwFileName[2002];
68    char structureFileName[2002];
69    char temp[200];
70    char readLine[500];
71 +  char xyzFile[200];
72 +  char shapeFile[200];
73 +  char shapeName[80];
74    FILE *vdwFile, *structureFile;
75    char* ffPath_env = "VDW_PATH";
76    char* ffPath;
# Line 64 | Line 78 | int main(int argc, char* argv[]){
78    char* foo;
79    char* myType;
80    char* vType;
81 +  char *token;
82 +  const char *file;
83 +  const char *period = ".";
84 +  int k;
85    int lineNum;
86    int nTokens;
87    int FF;
88 +  int bandwidth;
89 +  int gridwidth;
90    short int gotMatch;
91  
92    //parse the command line options
# Line 81 | Line 101 | int main(int argc, char* argv[]){
101      exit(1);
102    }
103  
104 +  file = fileName.c_str();
105 +  strcpy(xyzFile, file);
106 +  token = strtok(xyzFile, period);
107 +  strcpy(xyzFile, token);
108 +  strcpy(shapeFile, token);
109 +  strcpy(shapeName, token);
110 +  strcat(xyzFile, "Ref.xyz");
111 +  strcat(shapeFile, ".shape");
112 +  strcat(shapeName, "_RB_0");
113 +  ofstream xfiles(xyzFile);
114 +  
115 +  //the bandwidth has a default value (default=16), so it is always given
116 +  bandwidth = args_info.bandwidth_arg;
117 +  gridwidth = bandwidth*2;
118 +  
119    if (args_info.charmm_given) {
120      FF=CHARMM;
121      strcpy(vdwFileName, "charmm27.vdw");
# Line 106 | Line 141 | int main(int argc, char* argv[]){
141      strcpy(vdwFileName, "oplsaal.vdw");
142    }
143  
144 +
145 +  printf ("opening %s\n", vdwFileName);
146    vdwFile = fopen( vdwFileName, "r" );
147    
148    if( vdwFile == NULL ){
# Line 122 | Line 159 | int main(int argc, char* argv[]){
159      strcat( temp, vdwFileName );
160      strcpy( vdwFileName, temp );
161      
162 +    printf ("opening %s\n", vdwFileName);
163      vdwFile = fopen( vdwFileName, "r" );
164      
165      if( vdwFile == NULL ){
# Line 132 | Line 170 | int main(int argc, char* argv[]){
170                 vdwFileName );
171        exit(-1);
172      }
173 <    printf( "VDW file %s opened sucessfully.\n", vdwFileName );
174 <    lineNum = 0;
175 <    
173 >  }
174 >  printf( "VDW file %s opened sucessfully.\n", vdwFileName );
175 >  lineNum = 0;
176 >  
177 >  eof_test = fgets( readLine, sizeof(readLine), vdwFile );
178 >  lineNum++;
179 >  
180 >  if( eof_test == NULL ){
181 >    printf("Error in reading Atoms from force file at line %d.\n",
182 >           lineNum );
183 >    exit(-1);
184 >  }
185 >  
186 >  while( eof_test != NULL ){
187 >    // toss comment lines
188 >    if( readLine[0] != '!' && readLine[0] != '#' ){
189 >      
190 >      nTokens = count_tokens(readLine, " ,;\t");
191 >      if (nTokens < 4) {
192 >        printf("Not enough tokens on line %d.\n", lineNum);
193 >        exit(-1);
194 >      }
195 >      
196 >      at = new Atype();
197 >      foo = strtok(readLine, " ,;\t");
198 >      at->setType(foo);
199 >      foo = strtok(NULL, " ,;\t");      
200 >      mass = atof(foo);
201 >      at->setMass(mass);
202 >      foo = strtok(NULL, " ,;\t");
203 >      rpar = atof(foo);
204 >      at->setRpar(rpar);
205 >      foo = strtok(NULL, " ,;\t");
206 >      eps = atof(foo);  
207 >      at->setEps(eps);
208 >      vdwAtypes.push_back(at);        
209 >    }
210      eof_test = fgets( readLine, sizeof(readLine), vdwFile );
211      lineNum++;
212 <    
213 <    if( eof_test == NULL ){
214 <      printf("Error in reading Atoms from force file at line %d.\n",
215 <             lineNum );
144 <      exit(-1);
145 <    }
146 <    
147 <    while( eof_test != NULL ){
148 <      // toss comment lines
149 <      if( readLine[0] != '!' && readLine[0] != '#' ){
150 <      
151 <        nTokens = count_tokens(readLine, " ,;\t");
152 <        if (nTokens < 4) {
153 <          printf("Not enough tokens on line %d.\n", lineNum);
154 <          exit(-1);
155 <        }
156 <              
157 <        at = new Atype();
158 <        foo = strtok(readLine, " ,;\t");
159 <        at->setType(foo);
160 <        foo = strtok(NULL, " ,;\t");      
161 <        mass = atof(foo);
162 <        at->setMass(mass);
163 <        foo = strtok(NULL, " ,;\t");
164 <        rpar = atof(foo);
165 <        at->setRpar(rpar);
166 <        foo = strtok(NULL, " ,;\t");
167 <        eps = atof(foo);  
168 <        at->setEps(eps);
169 <        vdwAtypes.push_back(at);        
170 <      }
171 <      eof_test = fgets( readLine, sizeof(readLine), vdwFile );
172 <      lineNum++;
173 <    }          
174 <    
175 <    fclose( vdwFile );
176 <  }
212 >  }          
213 >  
214 >  fclose( vdwFile );
215 >  
216    printf("%d Atom Types read from VDW file\n", vdwAtypes.size());
217  
218    // Now, open and parse the input file!
# Line 222 | Line 261 | int main(int argc, char* argv[]){
261      }
262    }
263  
264 <  //GridBuilder gb = new GridBuilder();
265 <  //gb->findAxesAndOrigin(theAtoms);
266 <  //gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
267 <
264 >  rb = new RigidBody();
265 >  for( j = theAtoms.begin(); j !=  theAtoms.end(); ++j){
266 >    rb->addAtom(*j);
267 >  }
268  
269 +  rb->calcRefCoords();
270 +  
271 +  //print a reference coordinate xyz file
272 +  xfiles << rb->getNumAtoms() << "\n\n";
273 +  for (k=0; k<rb->getNumAtoms(); k++){
274 +    rb->getAtomRefCoor(xyz3, k);
275 +    xfiles << rb->getAtomBase(k) << "\t" <<
276 +              xyz3[0] << "\t" << xyz3[1] << "\t" <<
277 +              xyz3[2] << "\n";
278 +  }
279 +  
280 +  gb = new GridBuilder(rb, gridwidth);
281 +
282 +  cout << "Doing GridBuilder calculations...\n";
283 +  gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
284 +  
285 +  gb->printGridFiles();
286 +  
287 +  //load the grid element values to the main grid vectors
288 +  for (k=0; k<gridwidth*gridwidth; k++){
289 +    sigmaGrid.push_back(gb->passSig(k));
290 +    sGrid.push_back(gb->passS(k));
291 +    epsGrid.push_back(gb->passEps(k));
292 +  }
293 +
294 +  //build the spherical harmonic transform object
295 +  harmonize = new SphereHarm( bandwidth );
296 +  rb->getI(moments);
297 +  harmonize->printShapesFileStart(shapeFile, shapeName, rb->getMass(),
298 +                                  moments);
299 +  
300 +  //do the transforms and write to the shapes file
301 +  harmonize->doTransforms(sigmaGrid);
302 +  harmonize->printToShapesFile(shapeFile, 0);
303 +  harmonize->doTransforms(sGrid);
304 +  harmonize->printToShapesFile(shapeFile, 1);
305 +  harmonize->doTransforms(epsGrid);
306 +  harmonize->printToShapesFile(shapeFile, 2);
307 +  
308 +  //clean everything up
309 +  harmonize->~SphereHarm();  
310   }
311  
312   int count_tokens(char *line, char *delimiters) {

Diff Legend

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