trunk/SHAPES/shaper.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "shaperCmd.h"
#include "PDBReader.hpp"
#include "RigidBody.hpp"
#include "GridBuilder.hpp"
#include "SphereHarm.hpp"

#define MK_STR(s) # s
#define STR_DEFINE(t, s) t = MK_STR(s)
#define CHARMM 1
#define AMBER  2
#define LJ     3
#define GAFF   4
#define OPLS   5

int count_tokens(char *line, char *delimiters);
using namespace std;

class Atype{
public:
  Atype() {
    mass = 0.0;
    rpar = 0.0;
    eps = 0.0;
  }

  void setMass(double m) {mass = m;}
  void setRpar(double rp) {rpar = rp;}
  void setEps(double e) {eps = e;}
  void setType(char* t) {strcpy(type, t);}
  double getMass() {return mass;}
  double getRpar() {return rpar;}
  double getEps() {return eps;}
  char *getType() {return type;}
  
private:
  char type[100];
  double mass;
  double rpar;
  double eps;
};
  
int main(int argc, char* argv[]){

  gengetopt_args_info args_info;

  vector<Atype*>   vdwAtypes;
  vector<Atype*>::iterator i;
  Atype* at;
  RigidBody* rb;
  vector<VDWAtom*> theAtoms;
  vector<VDWAtom*>::iterator j;
  VDWAtom* atom;
  GridBuilder* gb;
  vector<double> sigmaGrid;
  vector<double> epsGrid;
  vector<double> sGrid;
  SphereHarm* harmonize;
  
  double mass, rpar, eps;
  double xyz3[3];
  double moments[3][3];
  double tolerance;
  string fileName;  
  char vdwFileName[2002];
  char structureFileName[2002];
  char temp[200];
  char readLine[500];
  char xyzFile[200];
  char shapeFile[200];
  char shapeName[80];
  FILE *vdwFile, *structureFile;
  char* ffPath_env = "VDW_PATH";
  char* ffPath;
  char* eof_test;
  char* foo;
  char* myType;
  char* vType;
  char *token;
  const char *file;
  const char *period = ".";
  int k;
  int lineNum;
  int nTokens;
  int FF;
  int bandwidth;
  int gridwidth;
  short int gotMatch;

  //parse the command line options
  if (cmdline_parser (argc, argv, &args_info) != 0)
    exit(1) ;

  if (args_info.input_given){
    fileName = args_info.input_arg;
  }
  else{
    std::cerr << "Does not have input file name" << endl;
    exit(1);
  }

  file = fileName.c_str();
  strcpy(xyzFile, file);
  token = strtok(xyzFile, period);
  strcpy(xyzFile, token);
  strcpy(shapeFile, token);
  strcpy(shapeName, token);
  strcat(xyzFile, "Ref.xyz");
  strcat(shapeFile, ".shape");
  ofstream xfiles(xyzFile);
  
  //the bandwidth has a default value (default=16), so it is always given
  bandwidth = args_info.bandwidth_arg;
  gridwidth = bandwidth*2;
 
  //the tolerance has a default value (default=0.01), so it is always given
  tolerance = args_info.tolerance_arg;

  if (args_info.charmm_given) {
    FF=CHARMM;
    strcpy(vdwFileName, "charmm27.vdw");
  }

  if (args_info.amber_given) {
    FF=AMBER;
    strcpy(vdwFileName, "amber99.vdw");
  }

  if (args_info.lj_given) {
    FF=LJ;
    strcpy(vdwFileName, "LJ.vdw");
  }

  if (args_info.gaff_given) {
    FF=GAFF;
    strcpy(vdwFileName, "gaff.vdw");
  }
  
  if (args_info.opls_given) {
    FF=OPLS;
    strcpy(vdwFileName, "oplsaal.vdw");
  }


  printf ("opening %s\n", vdwFileName);
  vdwFile = fopen( vdwFileName, "r" );
  
  if( vdwFile == NULL ){
       
    // next see if the force path enviorment variable is set
    
    ffPath = getenv( ffPath_env );
    if( ffPath == NULL ) {
      STR_DEFINE(ffPath, FRC_PATH );
    }
    
    strcpy( temp, ffPath );
    strcat( temp, "/" );
    strcat( temp, vdwFileName );
    strcpy( vdwFileName, temp );
    
    printf ("opening %s\n", vdwFileName);
    vdwFile = fopen( vdwFileName, "r" );
    
    if( vdwFile == NULL ){
      
      printf(  "Error opening the vanDerWaals parameter file: %s\n"
               "Have you tried setting the VDW_PARAM_DIR environment "
               "variable?\n",
               vdwFileName );
      exit(-1);
    }
  }
  printf( "VDW file %s opened sucessfully.\n", vdwFileName );
  lineNum = 0;
  
  eof_test = fgets( readLine, sizeof(readLine), vdwFile );
  lineNum++;
  
  if( eof_test == NULL ){
    printf("Error in reading Atoms from force file at line %d.\n",
           lineNum );
    exit(-1);
  }
  
  while( eof_test != NULL ){
    // toss comment lines
    if( readLine[0] != '!' && readLine[0] != '#' ){
      
      nTokens = count_tokens(readLine, " ,;\t");
      if (nTokens < 4) {
        printf("Not enough tokens on line %d.\n", lineNum);
        exit(-1);
      } 
      
      at = new Atype();
      foo = strtok(readLine, " ,;\t");
      at->setType(foo);
      foo = strtok(NULL, " ,;\t");       
      mass = atof(foo);
      at->setMass(mass);
      foo = strtok(NULL, " ,;\t");
      rpar = atof(foo);
      at->setRpar(rpar);
      foo = strtok(NULL, " ,;\t");
      eps = atof(foo);  
      at->setEps(eps);
      vdwAtypes.push_back(at);        
    }
    eof_test = fgets( readLine, sizeof(readLine), vdwFile );
    lineNum++;
  }          
  
  fclose( vdwFile );
  
  printf("%d Atom Types read from VDW file\n", vdwAtypes.size());

  // Now, open and parse the input file!

  strcpy(structureFileName, fileName.c_str() );

  PDBReader* PDBread = new PDBReader();
  PDBread->setPDBfile(structureFileName);
  theAtoms = PDBread->getAtomList();
  printf("Found %d atoms\n", theAtoms.size());

  for( j = theAtoms.begin(); j !=  theAtoms.end(); ++j){
    atom = *j;
    myType = atom->getType();
    gotMatch = 0;

    for( i = vdwAtypes.begin(); i !=  vdwAtypes.end(); ++i){
      at = *i;
      vType = at->getType();
      
      if (!strcasecmp(myType, vType)) {
        atom->setMass(at->getMass());
        atom->setRpar(at->getRpar());
        atom->setEps(at->getEps());
        gotMatch = 1;        
      }
    }
    if (!gotMatch) {      
      printf("No matches found for %s, ", myType);
      myType = atom->getBase();
      printf("trying with BaseType %s\n", myType);
      for( i = vdwAtypes.begin(); i !=  vdwAtypes.end(); ++i){
        at = *i;
        vType = at->getType();        
        if (!strcasecmp(myType, vType)) {
          atom->setMass(at->getMass());
          atom->setRpar(at->getRpar());
          atom->setEps(at->getEps());
          gotMatch = 1;
        }
      }
    }

    if (!gotMatch) {
      printf("No matches found with BaseType!\n");
    }
  }

  rb = new RigidBody();
  for( j = theAtoms.begin(); j !=  theAtoms.end(); ++j){
    rb->addAtom(*j);
  }

  rb->calcRefCoords();
  
  //print a reference coordinate xyz file
  xfiles << rb->getNumAtoms() << "\n\n";
  for (k=0; k<rb->getNumAtoms(); k++){
    rb->getAtomRefCoor(xyz3, k);
    xfiles << rb->getAtomBase(k) << "\t" <<
              xyz3[0] << "\t" << xyz3[1] << "\t" <<
              xyz3[2] << "\n";
  }
  
  gb = new GridBuilder(rb, gridwidth);

  cout << "Doing GridBuilder calculations...\n";
  gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
  
  gb->printGridFiles();
  
  //load the grid element values to the main grid vectors
  for (k=0; k<gridwidth*gridwidth; k++){
    sigmaGrid.push_back(gb->passSig(k));
    sGrid.push_back(gb->passS(k));
    epsGrid.push_back(gb->passEps(k));
  }

  //build the spherical harmonic transform object
  harmonize = new SphereHarm( bandwidth );
  rb->getI(moments);
  harmonize->printShapesFileStart(shapeFile, shapeName, rb->getMass(),
                                  moments);
  
  printf("Doing SHAPE calculations and outputting results...\n");
  //do the transforms and write to the shapes file
  harmonize->doTransforms(sigmaGrid);
  harmonize->printToShapesFile(shapeFile, 0, tolerance);
  harmonize->doTransforms(sGrid);
  harmonize->printToShapesFile(shapeFile, 1, tolerance);
  harmonize->doTransforms(epsGrid);
  harmonize->printToShapesFile(shapeFile, 2, tolerance);
  
  //clean everything up
  harmonize->~SphereHarm();  
  printf("Thank you for using SHAPES.  Enjoy your potential!\n");
}

int count_tokens(char *line, char *delimiters) {
/* PURPOSE: RETURN A COUNT OF THE NUMBER OF TOKENS ON THE LINE. */
 
  char *working_line;   /* WORKING COPY OF LINE. */
  int ntokens;          /* NUMBER OF TOKENS FOUND IN LINE. */
  char *strtok_ptr;     /* POINTER FOR STRTOK. */
   
  strtok_ptr= working_line= strdup(line);
   
  ntokens=0;
  while (strtok(strtok_ptr,delimiters)!=NULL)
    {
      ntokens++;
      strtok_ptr=NULL;
    }
   
  free(working_line);
  return(ntokens);
}
Revision:	1316
Committed:	Mon Jun 28 23:07:18 2004 UTC (20 years ago) by chrisfen
File size:	8254 byte(s)
Log Message:	Got rid of the _RB_0 in the shape name
#	User	Rev	Content
1	gezelter	1295	#include <iostream>
2			#include <fstream>
3			#include <string>
4			#include <vector>
5			#include "shaperCmd.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)
13			#define CHARMM 1
14			#define AMBER 2
15			#define LJ 3
16			#define GAFF 4
17			#define OPLS 5
18
19			int count_tokens(char line, char delimiters);
20			using namespace std;
21
22			class Atype{
23			public:
24			Atype() {
25			mass = 0.0;
26			rpar = 0.0;
27			eps = 0.0;
28			}
29
30			void setMass(double m) {mass = m;}
31			void setRpar(double rp) {rpar = rp;}
32			void setEps(double e) {eps = e;}
33			void setType(char* t) {strcpy(type, t);}
34			double getMass() {return mass;}
35			double getRpar() {return rpar;}
36			double getEps() {return eps;}
37			char *getType() {return type;}
38
39			private:
40			char type[100];
41			double mass;
42			double rpar;
43			double eps;
44			};
45
46			int main(int argc, char* argv[]){
47
48			gengetopt_args_info args_info;
49
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			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	chrisfen	1312	double tolerance;
67	gezelter	1295	string fileName;
68			char vdwFileName[2002];
69			char structureFileName[2002];
70			char temp[200];
71			char readLine[500];
72			char xyzFile[200];
73			char shapeFile[200];
74			char shapeName[80];
75			FILE vdwFile, structureFile;
76			char* ffPath_env = "VDW_PATH";
77			char* ffPath;
78			char* eof_test;
79			char* foo;
80			char* myType;
81			char* vType;
82			char *token;
83			const char *file;
84			const char *period = ".";
85			int k;
86			int lineNum;
87			int nTokens;
88			int FF;
89			int bandwidth;
90			int gridwidth;
91			short int gotMatch;
92
93			//parse the command line options
94			if (cmdline_parser (argc, argv, &args_info) != 0)
95			exit(1) ;
96
97			if (args_info.input_given){
98			fileName = args_info.input_arg;
99			}
100			else{
101			std::cerr << "Does not have input file name" << endl;
102			exit(1);
103			}
104
105			file = fileName.c_str();
106			strcpy(xyzFile, file);
107			token = strtok(xyzFile, period);
108			strcpy(xyzFile, token);
109			strcpy(shapeFile, token);
110			strcpy(shapeName, token);
111			strcat(xyzFile, "Ref.xyz");
112			strcat(shapeFile, ".shape");
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	chrisfen	1312
119			//the tolerance has a default value (default=0.01), so it is always given
120			tolerance = args_info.tolerance_arg;
121
122	gezelter	1295	if (args_info.charmm_given) {
123			FF=CHARMM;
124			strcpy(vdwFileName, "charmm27.vdw");
125			}
126
127			if (args_info.amber_given) {
128			FF=AMBER;
129			strcpy(vdwFileName, "amber99.vdw");
130			}
131
132			if (args_info.lj_given) {
133			FF=LJ;
134			strcpy(vdwFileName, "LJ.vdw");
135			}
136
137			if (args_info.gaff_given) {
138			FF=GAFF;
139			strcpy(vdwFileName, "gaff.vdw");
140			}
141
142			if (args_info.opls_given) {
143			FF=OPLS;
144			strcpy(vdwFileName, "oplsaal.vdw");
145			}
146
147
148			printf ("opening %s\n", vdwFileName);
149			vdwFile = fopen( vdwFileName, "r" );
150
151			if( vdwFile == NULL ){
152
153			// next see if the force path enviorment variable is set
154
155			ffPath = getenv( ffPath_env );
156			if( ffPath == NULL ) {
157			STR_DEFINE(ffPath, FRC_PATH );
158			}
159
160			strcpy( temp, ffPath );
161			strcat( temp, "/" );
162			strcat( temp, vdwFileName );
163			strcpy( vdwFileName, temp );
164
165			printf ("opening %s\n", vdwFileName);
166			vdwFile = fopen( vdwFileName, "r" );
167
168			if( vdwFile == NULL ){
169
170			printf( "Error opening the vanDerWaals parameter file: %s\n"
171			"Have you tried setting the VDW_PARAM_DIR environment "
172			"variable?\n",
173			vdwFileName );
174			exit(-1);
175			}
176			}
177			printf( "VDW file %s opened sucessfully.\n", vdwFileName );
178			lineNum = 0;
179
180			eof_test = fgets( readLine, sizeof(readLine), vdwFile );
181			lineNum++;
182
183			if( eof_test == NULL ){
184			printf("Error in reading Atoms from force file at line %d.\n",
185			lineNum );
186			exit(-1);
187			}
188
189			while( eof_test != NULL ){
190			// toss comment lines
191			if( readLine[0] != '!' && readLine[0] != '#' ){
192
193			nTokens = count_tokens(readLine, " ,;\t");
194			if (nTokens < 4) {
195			printf("Not enough tokens on line %d.\n", lineNum);
196			exit(-1);
197			}
198
199			at = new Atype();
200			foo = strtok(readLine, " ,;\t");
201			at->setType(foo);
202			foo = strtok(NULL, " ,;\t");
203			mass = atof(foo);
204			at->setMass(mass);
205			foo = strtok(NULL, " ,;\t");
206			rpar = atof(foo);
207			at->setRpar(rpar);
208			foo = strtok(NULL, " ,;\t");
209			eps = atof(foo);
210			at->setEps(eps);
211			vdwAtypes.push_back(at);
212			}
213			eof_test = fgets( readLine, sizeof(readLine), vdwFile );
214			lineNum++;
215			}
216
217			fclose( vdwFile );
218
219			printf("%d Atom Types read from VDW file\n", vdwAtypes.size());
220
221			// Now, open and parse the input file!
222
223			strcpy(structureFileName, fileName.c_str() );
224
225			PDBReader* PDBread = new PDBReader();
226			PDBread->setPDBfile(structureFileName);
227			theAtoms = PDBread->getAtomList();
228			printf("Found %d atoms\n", theAtoms.size());
229
230			for( j = theAtoms.begin(); j != theAtoms.end(); ++j){
231			atom = *j;
232			myType = atom->getType();
233			gotMatch = 0;
234
235			for( i = vdwAtypes.begin(); i != vdwAtypes.end(); ++i){
236			at = *i;
237			vType = at->getType();
238
239			if (!strcasecmp(myType, vType)) {
240			atom->setMass(at->getMass());
241			atom->setRpar(at->getRpar());
242			atom->setEps(at->getEps());
243			gotMatch = 1;
244			}
245			}
246			if (!gotMatch) {
247			printf("No matches found for %s, ", myType);
248			myType = atom->getBase();
249			printf("trying with BaseType %s\n", myType);
250			for( i = vdwAtypes.begin(); i != vdwAtypes.end(); ++i){
251			at = *i;
252			vType = at->getType();
253			if (!strcasecmp(myType, vType)) {
254			atom->setMass(at->getMass());
255			atom->setRpar(at->getRpar());
256			atom->setEps(at->getEps());
257			gotMatch = 1;
258			}
259			}
260			}
261
262			if (!gotMatch) {
263			printf("No matches found with BaseType!\n");
264			}
265			}
266
267			rb = new RigidBody();
268			for( j = theAtoms.begin(); j != theAtoms.end(); ++j){
269			rb->addAtom(*j);
270			}
271
272			rb->calcRefCoords();
273
274			//print a reference coordinate xyz file
275			xfiles << rb->getNumAtoms() << "\n\n";
276			for (k=0; k<rb->getNumAtoms(); k++){
277			rb->getAtomRefCoor(xyz3, k);
278			xfiles << rb->getAtomBase(k) << "\t" <<
279			xyz3[0] << "\t" << xyz3[1] << "\t" <<
280			xyz3[2] << "\n";
281			}
282
283			gb = new GridBuilder(rb, gridwidth);
284
285			cout << "Doing GridBuilder calculations...\n";
286			gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
287
288			gb->printGridFiles();
289
290			//load the grid element values to the main grid vectors
291			for (k=0; k<gridwidth*gridwidth; k++){
292			sigmaGrid.push_back(gb->passSig(k));
293			sGrid.push_back(gb->passS(k));
294			epsGrid.push_back(gb->passEps(k));
295			}
296
297			//build the spherical harmonic transform object
298			harmonize = new SphereHarm( bandwidth );
299			rb->getI(moments);
300			harmonize->printShapesFileStart(shapeFile, shapeName, rb->getMass(),
301			moments);
302
303			printf("Doing SHAPE calculations and outputting results...\n");
304			//do the transforms and write to the shapes file
305			harmonize->doTransforms(sigmaGrid);
306	chrisfen	1313	harmonize->printToShapesFile(shapeFile, 0, tolerance);
307	gezelter	1295	harmonize->doTransforms(sGrid);
308	chrisfen	1313	harmonize->printToShapesFile(shapeFile, 1, tolerance);
309	gezelter	1295	harmonize->doTransforms(epsGrid);
310	chrisfen	1313	harmonize->printToShapesFile(shapeFile, 2, tolerance);
311	gezelter	1295
312			//clean everything up
313			harmonize->~SphereHarm();
314			printf("Thank you for using SHAPES. Enjoy your potential!\n");
315			}
316
317			int count_tokens(char line, char delimiters) {
318			/* PURPOSE: RETURN A COUNT OF THE NUMBER OF TOKENS ON THE LINE. */
319
320			char working_line; / WORKING COPY OF LINE. */
321			int ntokens; /* NUMBER OF TOKENS FOUND IN LINE. */
322			char strtok_ptr; / POINTER FOR STRTOK. */
323
324			strtok_ptr= working_line= strdup(line);
325
326			ntokens=0;
327			while (strtok(strtok_ptr,delimiters)!=NULL)
328			{
329			ntokens++;
330			strtok_ptr=NULL;
331			}
332
333			free(working_line);
334			return(ntokens);
335			}