ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 33
Committed: Wed Jul 23 01:04:58 2008 UTC (11 years, 10 months ago) by tjod
File size: 14573 byte(s)
Log Message:
Don't write_sdf when missing parameters or other errors

Line File contents
1 #define EXTERN extern
2 #include "pcwin.h"
3 #include "pcmod.h"
4 #include "utility.h"
5 #include "energies.h"
6 #include "torsions.h"
7 #include "gmmx.h"
8 #include "pot.h"
9 #include "field.h"
10 #include "nonbond.h"
11
12 #include <time.h>
13 #include <errno.h>
14
15 struct t_logp {
16 float logp;
17 } logp_calc;
18 struct t_vibdata {
19 char ptgrp[4];
20 float mom_ix,mom_iy,mom_iz;
21 float etot,htot,stot,gtot,cptot;
22 } vibdata;
23
24 EXTERN struct t_optimize {
25 int param_avail, converge;
26 float initial_energy, final_energy, initial_heat, final_heat;
27 } optimize_data;
28
29 EXTERN struct t_files {
30 int nfiles, append, batch, icurrent, ibatno;
31 } files;
32 EXTERN struct t_dipolemom {
33 double total, xdipole, ydipole, zdipole;
34 } dipolemom;
35
36 EXTERN struct t_hfof {
37 float hfo, si;
38 int iunk;
39 } hfof;
40 EXTERN struct t_pcmfile {
41 char string[200];
42 int head;
43 char token[20];
44 int state;
45 unsigned int nocaps;
46 } pcmfile;
47 EXTERN struct t_minim_values {
48 int iprint, ndc, nconst;
49 float dielc;
50 } minim_values;
51 EXTERN struct t_minim_control {
52 int type, method, field, added_const;
53 char added_path[256],added_name[256];
54 } minim_control;
55 EXTERN struct t_rotbond {
56 int nbonds, bond[200][2], bond_res[200];
57 int incl_amide, incl_alkene;
58 } rotbond;
59
60 struct t_tmp {
61 int tmp_nrings, tmp_natoms[4], tmp_ratoms[30][4], tmp_clo_bond[4][4];
62 float tmp_clo_ang[2][4], tmp_clo_distance[4], tmp_ring_resolution[4];
63 } tmp;
64
65 int input_type;
66
67 int nerr;
68 int gmmx_abort;
69
70 int rd_sdf(FILE *);
71 int check_ring1(int);
72 int is_ring31(int);
73 int is_ring41(int);
74 int is_ring51(int);
75 int is_ring61(int);
76 int get_hybrid(int);
77 void set_atomtype(int,int,int,int,int,int);
78 void sort_file(char *, char *);
79 void gettoken(void);
80 int is_bond(int,int);
81 void read_atomtypes(int);
82 //double compute_rms(int, int **, double [150][3], double [150][3]);
83 //void fndumt( int, int **, double [150][3], double [150][3] );
84 //int fast_compare(double [150][3], double [150][3]);
85 //void slow_compare(int, int **, double [150][3], double [150][3]);
86 double compute_rms(int, int **, double **, double **);
87 void fndumt( int, int **, double **, double ** );
88 int fast_compare(double **, double **);
89 void slow_compare(int, int **, double **, double **);
90 void message_alert(char *, char *);
91 void gradient(void);
92 void eigen(int norder,double (*)[3],double *, double (*)[3]);
93 int check_structure(void);
94 int check_id(void);
95 double distance(int,int);
96 void rotb(int,int,int,int, float);
97 void make_bond(int , int , int );
98 void deletebond(int, int);
99 int setup_calculation(void);
100 void end_calculation(void);
101 void hster(int, int *);
102 int lepimer(int);
103 double dihdrl(int,int,int,int);
104 void mvatom(void);
105 void mvbond(void);
106 void minimize(void);
107 int findcd(int);
108 int findtor(int);
109 void InitialTransform(void);
110 int bad15(int);
111 int close_contact(int);
112 void gmmx2(void);
113 int check_stop(void);
114 void hbondreset(void);
115 void setup_output(char *);
116 void pimarkselection(void);
117 void mark_hbond(void);
118 void pireset(void);
119 void pcmfout(int);
120 void center(int);
121 void pcmfin(int,int);
122 void initialize(void);
123 void generate_bonds(void);
124 int bstereo(int);
125 void type(void);
126 void mmxsub(int);
127 void set_atom_color(void);
128 void eheat(void);
129 void charge_dipole(void);
130 void mmp22mod(int,int);
131 void check_numfile(int);
132 void dipole_dipole(void);
133 void set_active(void);
134 double energy(void);
135 void write_restart(void);
136 int read_restart(void);
137 void rdfile(int,int);
138 void montc(void);
139 void include_min(void);
140 void jacobi(int,int,double **,double *,double **,double *,double *);
141 void quatfit(int,int **,double **,double **);
142 FILE * fopen_path ( char * , char * , char * ) ;
143 void read_gmmxinp(char*,FILE*,int);
144 void run_gmmx(void);
145 int rmsck(int, int *);
146 void write_gmmx(void);
147 int read_sdf(int,int);
148 float ran1(int *);
149 void hdel(int);
150 void hadd(void);
151 void hcoord(void);
152 void inesc(char *);
153 void find_bonds(void);
154 int FetchRecord(FILE *, char *);
155 void deleteatom(int);
156 void write_sdf(int);
157 void zero_data(void);
158 void type_mmx(void);
159 void read_datafiles(char *);
160 void xlogp(float *);
161 void vibrate(void);
162
163 // =====================================
164 void read_gmmxinp(char *paramfile,
165 FILE *infile,
166 int flags)
167 {
168 char line[121];
169 int i,ncount,nret,j,k,jj,nc;
170 int icount;
171 int ngood, nbparam, nbatom,nbcontact;
172 int bcontact = FALSE;
173 int batom = FALSE;
174 float logp;
175 //char *std_file;
176 char bad_atom[80],bad_contact[80],bad_param[80];
177 clock_t start_time, end_time;
178 FILE *logfile;
179
180 // check filetype
181 ncount = 0;
182 input_type = 0;
183
184 /*
185 std_file = strdup(out_file);
186 strcpy(Savebox.fname, out_file);
187 */
188
189 // build filenames
190 strcpy(bad_atom,"_batom.sdf");
191 strcpy(bad_contact,"_bcontact.sdf");
192 strcpy(bad_param,"_bparam.sdf");
193
194 // open logfile
195 logfile = stderr;
196
197 files.nfiles = 0;
198 /*
199 // assume SDF input
200 while (FetchRecord(infile, line) )
201 {
202 if (strcmp(line,"$$$$") == 0)
203 files.nfiles++;
204 }
205 fclose(infile);
206 */
207 //
208 minim_control.field = MMFF94;
209 //strcpy(line,paramfile);
210 zero_data();
211 read_datafiles(paramfile);
212 gmmx.run = TRUE;
213 minim_control.method = 1; // just do first deriv minimization
214 icount = 0;
215 //
216 ngood = nbparam = nbatom = nbcontact = 0;
217 i = 0;
218 while (1)
219 // for (i=1; i <= files.nfiles; i++)
220 {
221 ++i;
222 initialize();
223 nret = rd_sdf(infile);
224 if (-1 == nret) return ;
225 if (FALSE == nret) {
226 fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
227 goto L_10;
228 natom = 0;
229 }
230
231 if (natom == 1)
232 {
233 fprintf(logfile,"Strnum: %d has only one atom\n",i);
234 goto L_10;
235 }
236 // if (i%20 == 0)
237 // printf("Strnum: %d\n",files.nfiles);
238 if (nret == TRUE)
239 {
240 nc = close_contact(1);
241 if (nc == FALSE)
242 {
243 fprintf(logfile,"Strnum: %d CID: %s fails on bad contact\n",files.nfiles,Struct_Title);
244 strcpy(Savebox.fname,bad_contact);
245 sprintf(Struct_Title,"Bad Contact %d\n",i);
246 files.append = TRUE;
247 if (bcontact == FALSE)
248 {
249 files.append = FALSE;
250 bcontact = TRUE;
251 }
252 nbcontact++;
253 //write_sdf(flags);
254 //strcpy(Savebox.fname,std_file);
255 goto L_10;
256 }
257 ncount = natom;
258 for (j=1; j<= ncount; j++)
259 {
260 jj = 0;
261 if (atom[j].atomnum == 1) // bridging hydrogens
262 {
263 if (atom[j].iat[0] != 0 && atom[j].iat[1] != 0)
264 goto L_10;
265 }
266 }
267 // TJO
268 if (flags & DO_ADDH) {
269 hdel(0);
270 hadd();
271 }
272 type();
273 // TJO
274 if (flags & DO_ADDH) {
275 hdel(0);
276 }
277 start_time = clock();
278 nret = setup_calculation();
279 if (nret == TRUE)
280 minimize();
281 end_calculation();
282
283 // minim_control.method = 3;
284 // TJO
285 if (flags & DO_ADDH) {
286 hdel(0);
287 hadd();
288 }
289 type_mmx();
290
291 nret = setup_calculation();
292 if (nret == TRUE)
293 {
294 minimize();
295 files.append = TRUE;
296 if ( icount == 0)
297 {
298 files.append = FALSE;
299 icount++;
300 }
301 // compute xlogp
302 if (energies.total < 1000.0)
303 {
304 if (flags & DO_XLOGP) {
305 if (VERBOSE) fprintf(stderr, " Evaluating XlogP\n");
306 xlogp(&logp);
307 logp_calc.logp = logp;
308 }
309 // compute dipole moment
310 dipolemom.xdipole = 0.0;
311 dipolemom.ydipole = 0.0;
312 dipolemom.zdipole = 0.0;
313 if (flags & DO_DIPOLE) {
314 if (VERBOSE) fprintf(stderr, " Evaluating dipole moment\n");
315 charge_dipole();
316 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
317 dipolemom.ydipole*dipolemom.ydipole +
318 dipolemom.zdipole*dipolemom.zdipole);
319 }
320 // compute vib
321 if (flags & DO_VIBRATION) {
322 if (VERBOSE) fprintf(stderr, " Evaluating vibrational parameters\n");
323 vibrate();
324 }
325 write_sdf(flags);
326 ++ngood;
327 ++files.nfiles;
328 }
329 } else
330 {
331 strcpy(Savebox.fname,bad_param);
332 files.append = TRUE;
333 //write_sdf(flags);
334 nbparam++;
335 //strcpy(Savebox.fname,std_file);
336 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
337 if (VERBOSE)
338 fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
339 }
340 end_calculation();
341 end_time = clock();
342 fprintf(logfile,"Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
343 natom,(end_time-start_time)/CLOCKS_PER_SEC);
344 if (VERBOSE) {
345 fprintf(stdout,
346 "Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
347 natom,(end_time-start_time)/CLOCKS_PER_SEC);
348 }
349 } else // read_sdf failed
350 {
351 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
352 strcpy(Savebox.fname,bad_atom);
353 files.append = TRUE;
354 if (batom == FALSE)
355 {
356 batom = TRUE;
357 files.append = FALSE;
358 }
359 //write_sdf(flags);
360 nbatom++;
361 //strcpy(Savebox.fname,std_file);
362 }
363 L_10:
364 continue;
365 }
366 fclose(infile);
367 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
368 if (VERBOSE) {
369 fprintf(stdout,
370 "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
371 ngood,nbparam,nbatom,nbcontact);
372 }
373
374 fclose(logfile);
375 }
376 // ================================================
377 // ==============================
378 void initialize_gmmx()
379 {
380 int i, j, k, nheav, nhyd,stmp;
381 long int mask;
382
383 mask = (1L << 0);
384 tmp.tmp_nrings = 0;
385 gmmxring.nrings = 0;
386 for (i=0; i < 4; i++)
387 {
388 tmp.tmp_natoms[i] = 0;
389 for (j=0; j < 4; j++)
390 {
391 gmmx_data.rng_size[j] = 0;
392 gmmxring.nring_atoms[i] = 0;
393 for (i=0; i < 30; i++)
394 {
395 gmmx_data.ring_atoms[i][j] = 0 ;
396 tmp.tmp_ratoms[i][j] = 0 ;
397 }
398 gmmx_data.clo_distance[j] = 0;
399 tmp.tmp_clo_distance[j] = 0 ;
400 gmmx_data.clo_ang[0][j] = 0;
401 gmmx_data.clo_ang[1][j] = 0;
402 tmp.tmp_clo_ang[0][j] = 0;
403 tmp.tmp_clo_ang[1][j] = 0;
404 gmmx_data.ring_resolution[j] = 0.0;
405 tmp.tmp_ring_resolution[j] = 0.0;
406 gmmx_data.clo_bond[0][j] = 0 ;
407 gmmx_data.clo_bond[1][j] = 0 ;
408 gmmx_data.clo_bond[2][j] = 0 ;
409 gmmx_data.clo_bond[3][j] = 0 ;
410 gmmx_data.clo_bond[4][j] = 0 ;
411 gmmx_data.clo_bond[5][j] = 0 ;
412 tmp.tmp_clo_bond[0][j] = 0 ;
413 tmp.tmp_clo_bond[1][j] = 0 ;
414 tmp.tmp_clo_bond[2][j] = 0 ;
415 tmp.tmp_clo_bond[3][j] = 0 ;
416 }
417 }
418 stmp = rand();
419 if(stmp >= 10000)
420 stmp = stmp/10;
421 if(stmp >= 10000)
422 stmp = stmp/10;
423 if(stmp >= 10000)
424 stmp = stmp/10;
425 if(stmp >= 10000)
426 stmp = stmp/10;
427 if(stmp >= 10000)
428 stmp = stmp/10;
429 if(stmp >= 10000)
430 stmp = stmp/10;
431 if(stmp >= 10000)
432 stmp = stmp/10;
433 if(stmp >= 10000)
434 stmp = stmp/10;
435 sprintf(gmmx_data.jobname,"gmx%d",stmp);
436 strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
437
438 gmmx_data.method = 3;
439 gmmx_data.iseed = 71277;
440 gmmx_data.its = 990;
441 gmmx_data.lnant = TRUE;
442 gmmx_data.hybrid = TRUE;
443 gmmx_data.ecut = TRUE;
444 gmmx_data.chig = TRUE;
445 gmmx_data.bad15 = TRUE;
446 gmmx_data.nopi1 = FALSE;
447 for (i=1; i <= natom; i++)
448 {
449 if (atom[i].flags & mask)
450 {
451 gmmx_data.nopi1 = TRUE;
452 break;
453 }
454 }
455 gmmx_data.hbond = TRUE;
456 gmmx_data.heat = FALSE;
457
458 gmmx_data.comp = 0;
459 gmmx_data.qpmr = FALSE;
460 gmmx_data.qdist = FALSE;
461 gmmx_data.qang = FALSE;
462 gmmx_data.qdihed = FALSE;
463 gmmx_data.npmr = 0;
464 gmmx_data.ndist = 0;
465 gmmx_data.nang = 0;
466 gmmx_data.ndihed = 0;
467
468 gmmx_data.nrings = 0;
469 gmmx_data.nbonds = 0;
470 gmmx_data.ewindow = 3.5;
471 gmmx_data.ewindow2 = 3.0;
472 gmmx_data.boltz = 300.0;
473 gmmx_data.ecutoff = 0.100;
474 gmmx_data.bad15_cutoff = 3.00;
475
476 gmmx_data.nsrms = 0;
477 gmmx_data.nsrmsa = 0;
478 gmmx_data.comp_method = 0;
479 gmmx_data.comp_method2 = 1;
480 gmmx_data.ermsa = 0.100;
481 gmmx_data.crmsa = 0.250;
482 gmmx_data.restart = FALSE;
483 gmmx_data.include_file = FALSE;
484
485 nheav = 0;
486 nhyd = 0;
487 for (i=1; i <= natom; i++)
488 {
489 if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
490 nheav++;
491 else
492 nhyd++;
493 }
494
495 gmmx_data.kstop = 5;
496 gmmx_data.kmin = 5*nheav;
497 gmmx_data.kdup = 5*nheav/2;
498 if (gmmx_data.kdup > 50)
499 gmmx_data.kdup = 50;
500 gmmx_data.max_search = 100000;
501
502 gmmx_data.nrings = 0;
503 for (i=0; i < 4; i++)
504 gmmx_data.nring_bonds[i] = 0;
505 for (i=0; i < 30; i++)
506 {
507 gmmxring.nring_atoms[i] = 0;
508 for (j=0; j < 30; j++)
509 gmmxring.ring_atoms[i][j]=0;
510 for (j=0; j < 4; j++)
511 for (k=0; k < 3; k++)
512 gmmx_data.ring_bond_data[j][i][k] = 0;
513 }
514
515 }
516 /* =============================================== */
517 // get hybridization of atom - tetrahedral = 1, planar = 2, linear = 3
518 //
519 int get_hybrid(int ia)
520 {
521 int itype;
522
523 itype = atom[ia].mmx_type;
524 if (itype == 2 || itype == 3 || itype == 7 || itype == 9 || itype == 25 ||
525 itype == 26 || itype == 29 || itype == 30 || itype == 37 || itype == 38 ||
526 itype == 40 || itype == 57 )
527 return 2;
528 else if (itype == 4 || itype == 10)
529 return 3;
530 else
531 return 1;
532 }
533 // ==================================
534 int check_ring1(int ia)
535 {
536 if (is_ring61(ia)) return TRUE;
537 if (is_ring51(ia)) return TRUE;
538 if (is_ring31(ia)) return TRUE;
539 if (is_ring41(ia)) return TRUE;
540 return FALSE;
541 }
542 // =======================================
543 int close_contact(int mode)
544 {
545 int i, j;
546 double dx,dy,dz;
547
548 for (i=1; i < natom; i++)
549 {
550 for (j=i+1; j <= natom; j++)
551 {
552 dx = fabs(atom[i].x - atom[j].x);
553 dy = fabs(atom[i].y - atom[j].y);
554 dz = fabs(atom[i].z - atom[j].z);
555 if ( (dx + dy + dz) < 0.1)
556 {
557 return FALSE;
558 }
559 }
560 }
561 return TRUE;
562 }