ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/gmmx_run.c
Revision: 25
Committed: Tue Jul 8 15:38:24 2008 UTC (13 years, 2 months ago) by tjod
File size: 14392 byte(s)
Log Message:
Incoporate the parameter files mmff94.prm and mmxconst.prm
into the mengine code by turning them into header files.
Each line becomes a string array element that is parsed in
place of lines in the .prm files.  The .prm files are no
longer needed, but are not deleted from the repository.

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