trunk/SHAPES/forcer.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "forcerCmd.h"
#include "PDBReader.hpp"
#include "RigidBody.hpp"
#include "GridBuilder.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;
  
  double mass, rpar, eps;
  string fileName;  
  char vdwFileName[2002];
  char structureFileName[2002];
  char temp[200];
  char readLine[500];
  FILE *vdwFile, *structureFile;
  char* ffPath_env = "VDW_PATH";
  char* ffPath;
  char* eof_test;
  char* foo;
  char* myType;
  char* vType;
  int lineNum;
  int nTokens;
  int FF;
  int bandwidth;
  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);
  }

  //the bandwidth has a default value (default=16), so it is always given
  bandwidth = args_info.bandwidth_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();
  
  gb = new GridBuilder(rb, bandwidth);

  cout << "Doing GridBuilder calculations...\n";
  gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
  
  //write out the grid files
  gb->printGridFiles();
}

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:	1281
Committed:	Mon Jun 21 13:38:55 2004 UTC (20 years, 2 months ago) by chrisfen
File size:	6350 byte(s)
Log Message:	Now calculates energies and builds the grid files. There are still bugs, but it compiles and runs...
#	User	Rev	Content
1	gezelter	1240	#include <iostream>
2			#include <fstream>
3			#include <string>
4			#include <vector>
5			#include "forcerCmd.h"
6			#include "PDBReader.hpp"
7	gezelter	1271	#include "RigidBody.hpp"
8	chrisfen	1279	#include "GridBuilder.hpp"
9	gezelter	1240
10			#define MK_STR(s) # s
11			#define STR_DEFINE(t, s) t = MK_STR(s)
12			#define CHARMM 1
13			#define AMBER 2
14			#define LJ 3
15			#define GAFF 4
16			#define OPLS 5
17
18			int count_tokens(char line, char delimiters);
19			using namespace std;
20
21			class Atype{
22			public:
23			Atype() {
24			mass = 0.0;
25			rpar = 0.0;
26			eps = 0.0;
27			}
28
29			void setMass(double m) {mass = m;}
30			void setRpar(double rp) {rpar = rp;}
31			void setEps(double e) {eps = e;}
32			void setType(char* t) {strcpy(type, t);}
33			double getMass() {return mass;}
34			double getRpar() {return rpar;}
35			double getEps() {return eps;}
36			char *getType() {return type;}
37
38			private:
39			char type[100];
40			double mass;
41			double rpar;
42			double eps;
43			};
44
45			int main(int argc, char* argv[]){
46
47			gengetopt_args_info args_info;
48
49			vector<Atype*> vdwAtypes;
50			vector<Atype*>::iterator i;
51			Atype* at;
52	gezelter	1271	RigidBody* rb;
53	gezelter	1240	vector<VDWAtom*> theAtoms;
54			vector<VDWAtom*>::iterator j;
55			VDWAtom* atom;
56	chrisfen	1279	GridBuilder* gb;
57			vector<double> sigmaGrid;
58			vector<double> epsGrid;
59			vector<double> sGrid;
60
61	gezelter	1240	double mass, rpar, eps;
62			string fileName;
63			char vdwFileName[2002];
64			char structureFileName[2002];
65			char temp[200];
66			char readLine[500];
67			FILE vdwFile, structureFile;
68	gezelter	1244	char* ffPath_env = "VDW_PATH";
69	gezelter	1240	char* ffPath;
70			char* eof_test;
71			char* foo;
72			char* myType;
73			char* vType;
74			int lineNum;
75			int nTokens;
76			int FF;
77	chrisfen	1275	int bandwidth;
78	gezelter	1240	short int gotMatch;
79
80			//parse the command line options
81			if (cmdline_parser (argc, argv, &args_info) != 0)
82			exit(1) ;
83
84			if (args_info.input_given){
85			fileName = args_info.input_arg;
86			}
87			else{
88			std::cerr << "Does not have input file name" << endl;
89			exit(1);
90			}
91
92	chrisfen	1281	//the bandwidth has a default value (default=16), so it is always given
93	chrisfen	1275	bandwidth = args_info.bandwidth_arg;
94
95	gezelter	1240	if (args_info.charmm_given) {
96			FF=CHARMM;
97			strcpy(vdwFileName, "charmm27.vdw");
98			}
99
100			if (args_info.amber_given) {
101			FF=AMBER;
102			strcpy(vdwFileName, "amber99.vdw");
103			}
104
105			if (args_info.lj_given) {
106			FF=LJ;
107			strcpy(vdwFileName, "LJ.vdw");
108			}
109
110			if (args_info.gaff_given) {
111			FF=GAFF;
112			strcpy(vdwFileName, "gaff.vdw");
113			}
114
115			if (args_info.opls_given) {
116			FF=OPLS;
117			strcpy(vdwFileName, "oplsaal.vdw");
118			}
119
120	gezelter	1274
121			printf ("opening %s\n", vdwFileName);
122	gezelter	1240	vdwFile = fopen( vdwFileName, "r" );
123
124			if( vdwFile == NULL ){
125
126			// next see if the force path enviorment variable is set
127
128			ffPath = getenv( ffPath_env );
129			if( ffPath == NULL ) {
130			STR_DEFINE(ffPath, FRC_PATH );
131			}
132
133			strcpy( temp, ffPath );
134			strcat( temp, "/" );
135			strcat( temp, vdwFileName );
136			strcpy( vdwFileName, temp );
137
138	gezelter	1274	printf ("opening %s\n", vdwFileName);
139	gezelter	1240	vdwFile = fopen( vdwFileName, "r" );
140
141			if( vdwFile == NULL ){
142
143			printf( "Error opening the vanDerWaals parameter file: %s\n"
144	gezelter	1244	"Have you tried setting the VDW_PARAM_DIR environment "
145	gezelter	1240	"variable?\n",
146			vdwFileName );
147			exit(-1);
148			}
149	gezelter	1274	}
150			printf( "VDW file %s opened sucessfully.\n", vdwFileName );
151			lineNum = 0;
152
153			eof_test = fgets( readLine, sizeof(readLine), vdwFile );
154			lineNum++;
155
156			if( eof_test == NULL ){
157			printf("Error in reading Atoms from force file at line %d.\n",
158			lineNum );
159			exit(-1);
160			}
161
162			while( eof_test != NULL ){
163			// toss comment lines
164			if( readLine[0] != '!' && readLine[0] != '#' ){
165
166			nTokens = count_tokens(readLine, " ,;\t");
167			if (nTokens < 4) {
168			printf("Not enough tokens on line %d.\n", lineNum);
169			exit(-1);
170			}
171
172			at = new Atype();
173			foo = strtok(readLine, " ,;\t");
174			at->setType(foo);
175			foo = strtok(NULL, " ,;\t");
176			mass = atof(foo);
177			at->setMass(mass);
178			foo = strtok(NULL, " ,;\t");
179			rpar = atof(foo);
180			at->setRpar(rpar);
181			foo = strtok(NULL, " ,;\t");
182			eps = atof(foo);
183			at->setEps(eps);
184			vdwAtypes.push_back(at);
185			}
186	gezelter	1240	eof_test = fgets( readLine, sizeof(readLine), vdwFile );
187			lineNum++;
188	gezelter	1274	}
189
190			fclose( vdwFile );
191
192	gezelter	1240	printf("%d Atom Types read from VDW file\n", vdwAtypes.size());
193
194			// Now, open and parse the input file!
195
196			strcpy(structureFileName, fileName.c_str() );
197
198			PDBReader* PDBread = new PDBReader();
199			PDBread->setPDBfile(structureFileName);
200			theAtoms = PDBread->getAtomList();
201			printf("Found %d atoms\n", theAtoms.size());
202
203			for( j = theAtoms.begin(); j != theAtoms.end(); ++j){
204			atom = *j;
205			myType = atom->getType();
206			gotMatch = 0;
207
208			for( i = vdwAtypes.begin(); i != vdwAtypes.end(); ++i){
209			at = *i;
210			vType = at->getType();
211
212			if (!strcasecmp(myType, vType)) {
213			atom->setMass(at->getMass());
214			atom->setRpar(at->getRpar());
215			atom->setEps(at->getEps());
216			gotMatch = 1;
217			}
218			}
219			if (!gotMatch) {
220			printf("No matches found for %s, ", myType);
221			myType = atom->getBase();
222			printf("trying with BaseType %s\n", myType);
223			for( i = vdwAtypes.begin(); i != vdwAtypes.end(); ++i){
224			at = *i;
225			vType = at->getType();
226			if (!strcasecmp(myType, vType)) {
227			atom->setMass(at->getMass());
228			atom->setRpar(at->getRpar());
229			atom->setEps(at->getEps());
230			gotMatch = 1;
231			}
232			}
233			}
234
235			if (!gotMatch) {
236			printf("No matches found with BaseType!\n");
237			}
238			}
239
240	gezelter	1271	rb = new RigidBody();
241			for( j = theAtoms.begin(); j != theAtoms.end(); ++j){
242	gezelter	1274	rb->addAtom(*j);
243	gezelter	1271	}
244	chrisfen	1281
245	chrisfen	1279	rb->calcRefCoords();
246
247			gb = new GridBuilder(rb, bandwidth);
248	chrisfen	1275
249	chrisfen	1281	cout << "Doing GridBuilder calculations...\n";
250	chrisfen	1279	gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
251	chrisfen	1281
252			//write out the grid files
253			gb->printGridFiles();
254	gezelter	1240	}
255
256			int count_tokens(char line, char delimiters) {
257			/* PURPOSE: RETURN A COUNT OF THE NUMBER OF TOKENS ON THE LINE. */
258
259			char working_line; / WORKING COPY OF LINE. */
260			int ntokens; /* NUMBER OF TOKENS FOUND IN LINE. */
261			char strtok_ptr; / POINTER FOR STRTOK. */
262
263			strtok_ptr= working_line= strdup(line);
264
265			ntokens=0;
266			while (strtok(strtok_ptr,delimiters)!=NULL)
267			{
268			ntokens++;
269			strtok_ptr=NULL;
270			}
271
272			free(working_line);
273			return(ntokens);
274			}