ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 75
Committed: Mon Dec 15 18:03:29 2008 UTC (12 years, 10 months ago) by gilbertke
File size: 9839 byte(s)
Log Message:
fixes in typeing to pass MMFF validation suite
Line File contents
1 /* NOTICE: this source code file has been modified for use with FreeMOL */
2 #define EXTERN extern
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "utility.h"
6 #include "energies.h"
7 #include "torsions.h"
8 #include "gmmx.h"
9 #include "pot.h"
10 #include "field.h"
11 #include "nonbond.h"
12 #include "dipmom.h"
13 #include "vibrate.h"
14
15 #include <errno.h>
16
17 struct t_logp {
18 float logp;
19 } logp_calc;
20
21 EXTERN struct t_files {
22 int nfiles, append, batch, icurrent, ibatno;
23 } files;
24
25 EXTERN struct t_pcmfile {
26 char string[200];
27 int head;
28 char token[20];
29 int state;
30 unsigned int nocaps;
31 } pcmfile;
32 EXTERN struct t_minim_values {
33 int iprint, ndc, nconst;
34 float dielc;
35 } minim_values;
36 EXTERN struct t_minim_control {
37 int type, method, field, added_const;
38 char added_path[256],added_name[256];
39 } minim_control;
40
41 struct t_tmp {
42 int tmp_nrings, tmp_natoms[4], tmp_ratoms[30][4], tmp_clo_bond[4][4];
43 float tmp_clo_ang[2][4], tmp_clo_distance[4], tmp_ring_resolution[4];
44 } tmp;
45
46 int input_type;
47
48 int nerr;
49 int gmmx_abort;
50
51 int rd_sdf(FILE *);
52
53 int setup_calculation(void);
54 void end_calculation(void);
55 void minimize(void);
56 static int close_contact(int);
57
58 void initialize(void);
59
60 void type(void);
61
62 void read_gmmxinp(char*,FILE*,int);
63 void hdel(int);
64 void hadd(void);
65 int FetchRecord(FILE *, char *);
66 void write_sdf(int);
67 void zero_data(void);
68 void type_mmx(void);
69 void read_datafiles(char *);
70 void xlogp(float *);
71
72 // =====================================
73 void read_gmmxinp(char *paramfile,
74 FILE *infile,
75 int flags)
76 {
77 char line[121];
78 int i,ncount,nret,j,k,jj,nc;
79 int icount;
80 int ngood, nbparam, nbatom,nbcontact;
81 int bcontact = FALSE;
82 int batom = FALSE;
83 float logp;
84 //char *std_file;
85 char bad_atom[80],bad_contact[80],bad_param[80];
86 FILE *logfile;
87
88 // check filetype
89 ncount = 0;
90 input_type = 0;
91
92 /*
93 std_file = strdup(out_file);
94 strcpy(Savebox.fname, out_file);
95 */
96
97 // build filenames
98 strcpy(bad_atom,"_batom.sdf");
99 strcpy(bad_contact,"_bcontact.sdf");
100 strcpy(bad_param,"_bparam.sdf");
101
102 // open logfile
103 logfile = pcmlogfile;
104
105 files.nfiles = 0;
106 /*
107 // assume SDF input
108 while (FetchRecord(infile, line) )
109 {
110 if (strcmp(line,"$$$$") == 0)
111 files.nfiles++;
112 }
113 fclose(infile);
114 */
115 //
116 minim_control.field = MMFF94;
117 //strcpy(line,paramfile);
118 zero_data();
119 read_datafiles(paramfile);
120 gmmx.run = TRUE;
121 minim_control.method = 1; // just do first deriv minimization
122 icount = 0;
123 //
124 ngood = nbparam = nbatom = nbcontact = 0;
125 i = 0;
126 while (1)
127 // for (i=1; i <= files.nfiles; i++)
128 {
129 ++i;
130 initialize();
131 nret = rd_sdf(infile);
132 if (-1 == nret) return ;
133 if (FALSE == nret)
134 {
135 fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
136 goto L_10;
137 natom = 0;
138 }
139
140 if (natom == 1)
141 {
142 fprintf(logfile,"Strnum: %d has only one atom\n",i);
143 goto L_10;
144 }
145 if (nret == TRUE)
146 {
147 // TJO
148 if (flags & DO_ADDH)
149 {
150 hdel(0);
151 hadd();
152 }
153 type();
154
155 // minim_values.iprint = TRUE;
156 nret = setup_calculation();
157 // for(i=1; i<= natom; i++)
158 //printf("Atom: %d %d %x\n",i,atom[i].type,atom[i].flags);
159 //if (minim_values.iprint == TRUE) exit(0);
160 if (nret == TRUE)
161 {
162 minimize();
163 // compute xlogp
164 if (energies.total < 1000.0)
165 {
166 if (flags & DO_XLOGP) {
167 if (VERBOSE) fprintf(pcmlogfile, " Evaluating XlogP\n");
168 xlogp(&logp);
169 logp_calc.logp = logp;
170 }
171 // compute dipole moment
172 dipolemom.xdipole = 0.0;
173 dipolemom.ydipole = 0.0;
174 dipolemom.zdipole = 0.0;
175 if (flags & DO_DIPOLE) {
176 if (VERBOSE) fprintf(pcmlogfile, " Evaluating dipole moment\n");
177 charge_dipole();
178 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
179 dipolemom.ydipole*dipolemom.ydipole +
180 dipolemom.zdipole*dipolemom.zdipole);
181 }
182 // compute vib
183 if (flags & DO_VIBRATION) {
184 if (VERBOSE) fprintf(pcmlogfile, " Evaluating vibrational parameters\n");
185 vibrate();
186 }
187 write_sdf(flags);
188 ++ngood;
189 ++files.nfiles;
190 }
191 } else
192 {
193 strcpy(Savebox.fname,bad_param);
194 files.append = TRUE;
195 //write_sdf(flags);
196 nbparam++;
197 //strcpy(Savebox.fname,std_file);
198 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
199 if (VERBOSE)
200 fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
201 }
202 end_calculation();
203 if (VERBOSE) {
204 fprintf(stdout,
205 "Strnum: %d %s natom %d \n",files.nfiles,Struct_Title, natom);
206 }
207 } else // read_sdf failed
208 {
209 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
210 strcpy(Savebox.fname,bad_atom);
211 files.append = TRUE;
212 if (batom == FALSE)
213 {
214 batom = TRUE;
215 files.append = FALSE;
216 }
217 //write_sdf(flags);
218 nbatom++;
219 //strcpy(Savebox.fname,std_file);
220 }
221 L_10:
222 continue;
223 }
224 // fclose(infile);
225 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
226 if (VERBOSE) {
227 fprintf(stdout,
228 "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
229 ngood,nbparam,nbatom,nbcontact);
230 }
231
232 // fclose(logfile);
233 }
234 // ================================================
235 // ==============================
236 void initialize_gmmx()
237 {
238 int i, j, k, nheav, nhyd,stmp;
239 long int mask;
240
241 mask = (1L << 0);
242 tmp.tmp_nrings = 0;
243 gmmxring.nrings = 0;
244 for (i=0; i < 4; i++)
245 {
246 tmp.tmp_natoms[i] = 0;
247 for (j=0; j < 4; j++)
248 {
249 gmmx_data.rng_size[j] = 0;
250 gmmxring.nring_atoms[i] = 0;
251 for (i=0; i < 30; i++)
252 {
253 gmmx_data.ring_atoms[i][j] = 0 ;
254 tmp.tmp_ratoms[i][j] = 0 ;
255 }
256 gmmx_data.clo_distance[j] = 0;
257 tmp.tmp_clo_distance[j] = 0 ;
258 gmmx_data.clo_ang[0][j] = 0;
259 gmmx_data.clo_ang[1][j] = 0;
260 tmp.tmp_clo_ang[0][j] = 0;
261 tmp.tmp_clo_ang[1][j] = 0;
262 gmmx_data.ring_resolution[j] = 0.0;
263 tmp.tmp_ring_resolution[j] = 0.0;
264 gmmx_data.clo_bond[0][j] = 0 ;
265 gmmx_data.clo_bond[1][j] = 0 ;
266 gmmx_data.clo_bond[2][j] = 0 ;
267 gmmx_data.clo_bond[3][j] = 0 ;
268 gmmx_data.clo_bond[4][j] = 0 ;
269 gmmx_data.clo_bond[5][j] = 0 ;
270 tmp.tmp_clo_bond[0][j] = 0 ;
271 tmp.tmp_clo_bond[1][j] = 0 ;
272 tmp.tmp_clo_bond[2][j] = 0 ;
273 tmp.tmp_clo_bond[3][j] = 0 ;
274 }
275 }
276 stmp = rand();
277 if(stmp >= 10000)
278 stmp = stmp/10;
279 if(stmp >= 10000)
280 stmp = stmp/10;
281 if(stmp >= 10000)
282 stmp = stmp/10;
283 if(stmp >= 10000)
284 stmp = stmp/10;
285 if(stmp >= 10000)
286 stmp = stmp/10;
287 if(stmp >= 10000)
288 stmp = stmp/10;
289 if(stmp >= 10000)
290 stmp = stmp/10;
291 if(stmp >= 10000)
292 stmp = stmp/10;
293 sprintf(gmmx_data.jobname,"gmx%d",stmp);
294 strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
295
296 gmmx_data.method = 3;
297 gmmx_data.iseed = 71277;
298 gmmx_data.its = 990;
299 gmmx_data.lnant = TRUE;
300 gmmx_data.hybrid = TRUE;
301 gmmx_data.ecut = TRUE;
302 gmmx_data.chig = TRUE;
303 gmmx_data.bad15 = TRUE;
304 gmmx_data.nopi1 = FALSE;
305 for (i=1; i <= natom; i++)
306 {
307 if (atom[i].flags & mask)
308 {
309 gmmx_data.nopi1 = TRUE;
310 break;
311 }
312 }
313 gmmx_data.hbond = TRUE;
314 gmmx_data.heat = FALSE;
315
316 gmmx_data.comp = 0;
317 gmmx_data.qpmr = FALSE;
318 gmmx_data.qdist = FALSE;
319 gmmx_data.qang = FALSE;
320 gmmx_data.qdihed = FALSE;
321 gmmx_data.npmr = 0;
322 gmmx_data.ndist = 0;
323 gmmx_data.nang = 0;
324 gmmx_data.ndihed = 0;
325
326 gmmx_data.nrings = 0;
327 gmmx_data.nbonds = 0;
328 gmmx_data.ewindow = 3.5;
329 gmmx_data.ewindow2 = 3.0;
330 gmmx_data.boltz = 300.0;
331 gmmx_data.ecutoff = 0.100;
332 gmmx_data.bad15_cutoff = 3.00;
333
334 gmmx_data.nsrms = 0;
335 gmmx_data.nsrmsa = 0;
336 gmmx_data.comp_method = 0;
337 gmmx_data.comp_method2 = 1;
338 gmmx_data.ermsa = 0.100;
339 gmmx_data.crmsa = 0.250;
340 gmmx_data.restart = FALSE;
341 gmmx_data.include_file = FALSE;
342
343 nheav = 0;
344 nhyd = 0;
345 for (i=1; i <= natom; i++)
346 {
347 if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
348 nheav++;
349 else
350 nhyd++;
351 }
352
353 gmmx_data.kstop = 5;
354 gmmx_data.kmin = 5*nheav;
355 gmmx_data.kdup = 5*nheav/2;
356 if (gmmx_data.kdup > 50)
357 gmmx_data.kdup = 50;
358 gmmx_data.max_search = 100000;
359
360 gmmx_data.nrings = 0;
361 for (i=0; i < 4; i++)
362 gmmx_data.nring_bonds[i] = 0;
363 for (i=0; i < 30; i++)
364 {
365 gmmxring.nring_atoms[i] = 0;
366 for (j=0; j < 30; j++)
367 gmmxring.ring_atoms[i][j]=0;
368 for (j=0; j < 4; j++)
369 for (k=0; k < 3; k++)
370 gmmx_data.ring_bond_data[j][i][k] = 0;
371 }
372
373 }
374 // =======================================
375 int close_contact(int mode)
376 {
377 int i, j;
378 double dx,dy,dz;
379
380 for (i=1; i < natom; i++)
381 {
382 for (j=i+1; j <= natom; j++)
383 {
384 dx = fabs(atom[i].x - atom[j].x);
385 dy = fabs(atom[i].y - atom[j].y);
386 dz = fabs(atom[i].z - atom[j].z);
387 if ( (dx + dy + dz) < 0.1)
388 {
389 return FALSE;
390 }
391 }
392 }
393 return TRUE;
394 }