ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 48
Committed: Thu Jul 31 16:34:52 2008 UTC (11 years, 8 months ago) by tjod
File size: 14648 byte(s)
Log Message:
Add modification notice to files modified for Freemol
per mengine/smi23d license requirements

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