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, 11 months ago) by tjod
File size: 14887 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #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     }