ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/gmmx_run.c
Revision: 16
Committed: Mon Jun 23 21:29:26 2008 UTC (11 years, 9 months ago) by tjod
File size: 14239 byte(s)
Log Message:
Add back the command line options -d -i -x for dipole moment, vibrational data and logp.  Have error.log directed to stdout.

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(line);
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 while (1)
218 // for (i=1; i <= files.nfiles; i++)
219 {
220 initialize();
221 nret = rd_sdf(infile);
222 if (-1 == nret) return ;
223
224 if (natom == 1)
225 {
226 fprintf(logfile,"Strnum: %d has only one atom\n",i);
227 goto L_10;
228 }
229 // if (i%20 == 0)
230 // printf("Strnum: %d\n",files.nfiles);
231 if (nret == TRUE)
232 {
233 nc = close_contact(1);
234 if (nc == FALSE)
235 {
236 fprintf(logfile,"Strnum: %d CID: %s fails on bad contact\n",files.nfiles,Struct_Title);
237 strcpy(Savebox.fname,bad_contact);
238 sprintf(Struct_Title,"Bad Contact %d\n",i);
239 files.append = TRUE;
240 if (bcontact == FALSE)
241 {
242 files.append = FALSE;
243 bcontact = TRUE;
244 }
245 nbcontact++;
246 write_sdf(flags);
247 //strcpy(Savebox.fname,std_file);
248 goto L_10;
249 }
250 ncount = natom;
251 for (j=1; j<= ncount; j++)
252 {
253 jj = 0;
254 if (atom[j].atomnum == 1) // bridging hydrogens
255 {
256 if (atom[j].iat[0] != 0 && atom[j].iat[1] != 0)
257 goto L_10;
258 }
259 }
260 hdel(0);
261 hadd();
262 type();
263 hdel(0);
264 start_time = clock();
265 nret = setup_calculation();
266 if (nret == TRUE)
267 minimize();
268 end_calculation();
269
270 // minim_control.method = 3;
271 hdel(0);
272 hadd();
273 type_mmx();
274
275 nret = setup_calculation();
276 if (nret == TRUE)
277 {
278 minimize();
279 files.append = TRUE;
280 if ( icount == 0)
281 {
282 files.append = FALSE;
283 icount++;
284 }
285 // compute xlogp
286 if (energies.total < 1000.0)
287 {
288 if (flags & DO_XLOGP) {
289 if (VERBOSE) fprintf(stderr, " Evaluating XlogP\n");
290 xlogp(&logp);
291 logp_calc.logp = logp;
292 }
293 // compute dipole moment
294 dipolemom.xdipole = 0.0;
295 dipolemom.ydipole = 0.0;
296 dipolemom.zdipole = 0.0;
297 if (flags & DO_DIPOLE) {
298 if (VERBOSE) fprintf(stderr, " Evaluating dipole moment\n");
299 charge_dipole();
300 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
301 dipolemom.ydipole*dipolemom.ydipole +
302 dipolemom.zdipole*dipolemom.zdipole);
303 }
304 // compute vib
305 if (flags & DO_VIBRATION) {
306 if (VERBOSE) fprintf(stderr, " Evaluating vibrational parameters\n");
307 vibrate();
308 }
309 write_sdf(flags);
310 ++ngood;
311 ++files.nfiles;
312 }
313 } else
314 {
315 strcpy(Savebox.fname,bad_param);
316 files.append = TRUE;
317 write_sdf(flags);
318 nbparam++;
319 //strcpy(Savebox.fname,std_file);
320 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
321 if (VERBOSE)
322 fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
323 }
324 end_calculation();
325 end_time = clock();
326 fprintf(logfile,"Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
327 natom,(end_time-start_time)/CLOCKS_PER_SEC);
328 if (VERBOSE) {
329 fprintf(stdout,
330 "Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
331 natom,(end_time-start_time)/CLOCKS_PER_SEC);
332 }
333 } else // read_sdf failed
334 {
335 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
336 strcpy(Savebox.fname,bad_atom);
337 files.append = TRUE;
338 if (batom == FALSE)
339 {
340 batom = TRUE;
341 files.append = FALSE;
342 }
343 write_sdf(flags);
344 nbatom++;
345 //strcpy(Savebox.fname,std_file);
346 }
347 L_10:
348 continue;
349 }
350 fclose(infile);
351 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
352 if (VERBOSE) {
353 fprintf(stdout,
354 "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
355 ngood,nbparam,nbatom,nbcontact);
356 }
357
358 fclose(logfile);
359 }
360 // ================================================
361 // ==============================
362 void initialize_gmmx()
363 {
364 int i, j, k, nheav, nhyd,stmp;
365 long int mask;
366
367 mask = (1L << 0);
368 tmp.tmp_nrings = 0;
369 gmmxring.nrings = 0;
370 for (i=0; i < 4; i++)
371 {
372 tmp.tmp_natoms[i] = 0;
373 for (j=0; j < 4; j++)
374 {
375 gmmx_data.rng_size[j] = 0;
376 gmmxring.nring_atoms[i] = 0;
377 for (i=0; i < 30; i++)
378 {
379 gmmx_data.ring_atoms[i][j] = 0 ;
380 tmp.tmp_ratoms[i][j] = 0 ;
381 }
382 gmmx_data.clo_distance[j] = 0;
383 tmp.tmp_clo_distance[j] = 0 ;
384 gmmx_data.clo_ang[0][j] = 0;
385 gmmx_data.clo_ang[1][j] = 0;
386 tmp.tmp_clo_ang[0][j] = 0;
387 tmp.tmp_clo_ang[1][j] = 0;
388 gmmx_data.ring_resolution[j] = 0.0;
389 tmp.tmp_ring_resolution[j] = 0.0;
390 gmmx_data.clo_bond[0][j] = 0 ;
391 gmmx_data.clo_bond[1][j] = 0 ;
392 gmmx_data.clo_bond[2][j] = 0 ;
393 gmmx_data.clo_bond[3][j] = 0 ;
394 gmmx_data.clo_bond[4][j] = 0 ;
395 gmmx_data.clo_bond[5][j] = 0 ;
396 tmp.tmp_clo_bond[0][j] = 0 ;
397 tmp.tmp_clo_bond[1][j] = 0 ;
398 tmp.tmp_clo_bond[2][j] = 0 ;
399 tmp.tmp_clo_bond[3][j] = 0 ;
400 }
401 }
402 stmp = rand();
403 if(stmp >= 10000)
404 stmp = stmp/10;
405 if(stmp >= 10000)
406 stmp = stmp/10;
407 if(stmp >= 10000)
408 stmp = stmp/10;
409 if(stmp >= 10000)
410 stmp = stmp/10;
411 if(stmp >= 10000)
412 stmp = stmp/10;
413 if(stmp >= 10000)
414 stmp = stmp/10;
415 if(stmp >= 10000)
416 stmp = stmp/10;
417 if(stmp >= 10000)
418 stmp = stmp/10;
419 sprintf(gmmx_data.jobname,"gmx%d",stmp);
420 strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
421
422 gmmx_data.method = 3;
423 gmmx_data.iseed = 71277;
424 gmmx_data.its = 990;
425 gmmx_data.lnant = TRUE;
426 gmmx_data.hybrid = TRUE;
427 gmmx_data.ecut = TRUE;
428 gmmx_data.chig = TRUE;
429 gmmx_data.bad15 = TRUE;
430 gmmx_data.nopi1 = FALSE;
431 for (i=1; i <= natom; i++)
432 {
433 if (atom[i].flags & mask)
434 {
435 gmmx_data.nopi1 = TRUE;
436 break;
437 }
438 }
439 gmmx_data.hbond = TRUE;
440 gmmx_data.heat = FALSE;
441
442 gmmx_data.comp = 0;
443 gmmx_data.qpmr = FALSE;
444 gmmx_data.qdist = FALSE;
445 gmmx_data.qang = FALSE;
446 gmmx_data.qdihed = FALSE;
447 gmmx_data.npmr = 0;
448 gmmx_data.ndist = 0;
449 gmmx_data.nang = 0;
450 gmmx_data.ndihed = 0;
451
452 gmmx_data.nrings = 0;
453 gmmx_data.nbonds = 0;
454 gmmx_data.ewindow = 3.5;
455 gmmx_data.ewindow2 = 3.0;
456 gmmx_data.boltz = 300.0;
457 gmmx_data.ecutoff = 0.100;
458 gmmx_data.bad15_cutoff = 3.00;
459
460 gmmx_data.nsrms = 0;
461 gmmx_data.nsrmsa = 0;
462 gmmx_data.comp_method = 0;
463 gmmx_data.comp_method2 = 1;
464 gmmx_data.ermsa = 0.100;
465 gmmx_data.crmsa = 0.250;
466 gmmx_data.restart = FALSE;
467 gmmx_data.include_file = FALSE;
468
469 nheav = 0;
470 nhyd = 0;
471 for (i=1; i <= natom; i++)
472 {
473 if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
474 nheav++;
475 else
476 nhyd++;
477 }
478
479 gmmx_data.kstop = 5;
480 gmmx_data.kmin = 5*nheav;
481 gmmx_data.kdup = 5*nheav/2;
482 if (gmmx_data.kdup > 50)
483 gmmx_data.kdup = 50;
484 gmmx_data.max_search = 100000;
485
486 gmmx_data.nrings = 0;
487 for (i=0; i < 4; i++)
488 gmmx_data.nring_bonds[i] = 0;
489 for (i=0; i < 30; i++)
490 {
491 gmmxring.nring_atoms[i] = 0;
492 for (j=0; j < 30; j++)
493 gmmxring.ring_atoms[i][j]=0;
494 for (j=0; j < 4; j++)
495 for (k=0; k < 3; k++)
496 gmmx_data.ring_bond_data[j][i][k] = 0;
497 }
498
499 }
500 /* =============================================== */
501 // get hybridization of atom - tetrahedral = 1, planar = 2, linear = 3
502 //
503 int get_hybrid(int ia)
504 {
505 int itype;
506
507 itype = atom[ia].mmx_type;
508 if (itype == 2 || itype == 3 || itype == 7 || itype == 9 || itype == 25 ||
509 itype == 26 || itype == 29 || itype == 30 || itype == 37 || itype == 38 ||
510 itype == 40 || itype == 57 )
511 return 2;
512 else if (itype == 4 || itype == 10)
513 return 3;
514 else
515 return 1;
516 }
517 // ==================================
518 int check_ring1(int ia)
519 {
520 if (is_ring61(ia)) return TRUE;
521 if (is_ring51(ia)) return TRUE;
522 if (is_ring31(ia)) return TRUE;
523 if (is_ring41(ia)) return TRUE;
524 return FALSE;
525 }
526 // =======================================
527 int close_contact(int mode)
528 {
529 int i, j;
530 double dx,dy,dz;
531
532 for (i=1; i < natom; i++)
533 {
534 for (j=i+1; j <= natom; j++)
535 {
536 dx = fabs(atom[i].x - atom[j].x);
537 dy = fabs(atom[i].y - atom[j].y);
538 dz = fabs(atom[i].z - atom[j].z);
539 if ( (dx + dy + dz) < 0.1)
540 {
541 return FALSE;
542 }
543 }
544 }
545 return TRUE;
546 }