trunk/SHAPES/forcer.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "forcerCmd.h"
#include "PDBReader.hpp"
#include "RigidBody.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;

  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;
  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);
  }

  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");
  }

  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 );
    
    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);
  }

  //GridBuilder gb = new GridBuilder();
  //gb->findAxesAndOrigin(theAtoms);
  //gb->launchProbe(FF, sigmaGrid, sGrid, epsGrid);
 

}

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