ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/gmmx_run.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 14887 byte(s)
Log Message:
test

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