ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 67
Committed: Fri Dec 5 14:18:02 2008 UTC (11 years, 9 months ago) by gilbertke
File size: 9849 byte(s)
Log Message:
added aromatic bonds in read_sdf removed TJO bailout
Line User Rev File contents
1 tjod 48 /* NOTICE: this source code file has been modified for use with FreeMOL */
2 tjod 3 #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 gilbertke 63 #include "dipmom.h"
13     #include "vibrate.h"
14 tjod 3
15     #include <errno.h>
16    
17     struct t_logp {
18     float logp;
19     } logp_calc;
20    
21     EXTERN struct t_files {
22     int nfiles, append, batch, icurrent, ibatno;
23     } files;
24    
25     EXTERN struct t_pcmfile {
26     char string[200];
27     int head;
28     char token[20];
29     int state;
30     unsigned int nocaps;
31     } pcmfile;
32     EXTERN struct t_minim_values {
33     int iprint, ndc, nconst;
34     float dielc;
35     } minim_values;
36     EXTERN struct t_minim_control {
37     int type, method, field, added_const;
38     char added_path[256],added_name[256];
39     } minim_control;
40    
41     struct t_tmp {
42     int tmp_nrings, tmp_natoms[4], tmp_ratoms[30][4], tmp_clo_bond[4][4];
43     float tmp_clo_ang[2][4], tmp_clo_distance[4], tmp_ring_resolution[4];
44     } tmp;
45    
46     int input_type;
47    
48     int nerr;
49     int gmmx_abort;
50    
51     int rd_sdf(FILE *);
52 gilbertke 63
53 tjod 3 int setup_calculation(void);
54     void end_calculation(void);
55     void minimize(void);
56 gilbertke 63 static int close_contact(int);
57    
58 tjod 3 void initialize(void);
59 gilbertke 63
60 tjod 3 void type(void);
61 gilbertke 63
62 tjod 16 void read_gmmxinp(char*,FILE*,int);
63 tjod 3 void hdel(int);
64     void hadd(void);
65     int FetchRecord(FILE *, char *);
66     void write_sdf(int);
67     void zero_data(void);
68     void type_mmx(void);
69     void read_datafiles(char *);
70     void xlogp(float *);
71    
72     // =====================================
73 tjod 16 void read_gmmxinp(char *paramfile,
74 tjod 15 FILE *infile,
75 tjod 3 int flags)
76     {
77     char line[121];
78     int i,ncount,nret,j,k,jj,nc;
79     int icount;
80     int ngood, nbparam, nbatom,nbcontact;
81     int bcontact = FALSE;
82     int batom = FALSE;
83     float logp;
84 tjod 15 //char *std_file;
85 tjod 3 char bad_atom[80],bad_contact[80],bad_param[80];
86 tjod 15 FILE *logfile;
87 tjod 3
88     // check filetype
89     ncount = 0;
90     input_type = 0;
91    
92 tjod 15 /*
93 tjod 3 std_file = strdup(out_file);
94     strcpy(Savebox.fname, out_file);
95 tjod 15 */
96 tjod 3
97     // build filenames
98 tjod 15 strcpy(bad_atom,"_batom.sdf");
99     strcpy(bad_contact,"_bcontact.sdf");
100     strcpy(bad_param,"_bparam.sdf");
101 tjod 3
102     // open logfile
103 wdelano 58 logfile = pcmlogfile;
104 tjod 3
105     files.nfiles = 0;
106 tjod 15 /*
107 tjod 3 // assume SDF input
108     while (FetchRecord(infile, line) )
109     {
110     if (strcmp(line,"$$$$") == 0)
111     files.nfiles++;
112     }
113     fclose(infile);
114 tjod 15 */
115 tjod 3 //
116     minim_control.field = MMFF94;
117 tjod 25 //strcpy(line,paramfile);
118 tjod 3 zero_data();
119 tjod 25 read_datafiles(paramfile);
120 tjod 3 gmmx.run = TRUE;
121     minim_control.method = 1; // just do first deriv minimization
122     icount = 0;
123     //
124     ngood = nbparam = nbatom = nbcontact = 0;
125 tjod 31 i = 0;
126 tjod 15 while (1)
127     // for (i=1; i <= files.nfiles; i++)
128 tjod 3 {
129 tjod 31 ++i;
130 tjod 3 initialize();
131     nret = rd_sdf(infile);
132 tjod 15 if (-1 == nret) return ;
133 gilbertke 64 if (FALSE == nret)
134     {
135 tjod 31 fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
136     goto L_10;
137     natom = 0;
138 gilbertke 64 }
139 tjod 3
140     if (natom == 1)
141     {
142     fprintf(logfile,"Strnum: %d has only one atom\n",i);
143     goto L_10;
144     }
145     if (nret == TRUE)
146     {
147 tjod 20 // TJO
148 gilbertke 64 if (flags & DO_ADDH)
149     {
150 tjod 20 hdel(0);
151     hadd();
152     }
153 gilbertke 65 type();
154 tjod 3
155 gilbertke 67 // minim_values.iprint = TRUE;
156 tjod 3 nret = setup_calculation();
157 gilbertke 67 // if (minim_values.iprint == TRUE) exit(0);
158 tjod 3 if (nret == TRUE)
159     {
160     minimize();
161     files.append = TRUE;
162     if ( icount == 0)
163     {
164     files.append = FALSE;
165     icount++;
166     }
167     // compute xlogp
168     if (energies.total < 1000.0)
169     {
170     if (flags & DO_XLOGP) {
171 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating XlogP\n");
172 tjod 3 xlogp(&logp);
173     logp_calc.logp = logp;
174     }
175     // compute dipole moment
176     dipolemom.xdipole = 0.0;
177     dipolemom.ydipole = 0.0;
178     dipolemom.zdipole = 0.0;
179     if (flags & DO_DIPOLE) {
180 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating dipole moment\n");
181 tjod 3 charge_dipole();
182     dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
183     dipolemom.ydipole*dipolemom.ydipole +
184     dipolemom.zdipole*dipolemom.zdipole);
185     }
186     // compute vib
187     if (flags & DO_VIBRATION) {
188 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating vibrational parameters\n");
189 tjod 3 vibrate();
190     }
191     write_sdf(flags);
192 tjod 15 ++ngood;
193     ++files.nfiles;
194 tjod 3 }
195     } else
196     {
197     strcpy(Savebox.fname,bad_param);
198     files.append = TRUE;
199 tjod 33 //write_sdf(flags);
200 tjod 3 nbparam++;
201 tjod 15 //strcpy(Savebox.fname,std_file);
202     fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
203 tjod 3 if (VERBOSE)
204     fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
205     }
206     end_calculation();
207     if (VERBOSE) {
208     fprintf(stdout,
209 gilbertke 63 "Strnum: %d %s natom %d \n",files.nfiles,Struct_Title, natom);
210 tjod 3 }
211     } else // read_sdf failed
212     {
213 tjod 15 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
214 tjod 3 strcpy(Savebox.fname,bad_atom);
215     files.append = TRUE;
216     if (batom == FALSE)
217     {
218     batom = TRUE;
219     files.append = FALSE;
220     }
221 tjod 33 //write_sdf(flags);
222 tjod 3 nbatom++;
223 tjod 15 //strcpy(Savebox.fname,std_file);
224 tjod 3 }
225     L_10:
226     continue;
227     }
228 wdelano 58 // fclose(infile);
229 tjod 3 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
230     if (VERBOSE) {
231     fprintf(stdout,
232     "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
233     ngood,nbparam,nbatom,nbcontact);
234     }
235    
236 wdelano 58 // fclose(logfile);
237 tjod 3 }
238     // ================================================
239     // ==============================
240     void initialize_gmmx()
241     {
242     int i, j, k, nheav, nhyd,stmp;
243     long int mask;
244    
245     mask = (1L << 0);
246     tmp.tmp_nrings = 0;
247     gmmxring.nrings = 0;
248     for (i=0; i < 4; i++)
249     {
250     tmp.tmp_natoms[i] = 0;
251     for (j=0; j < 4; j++)
252     {
253     gmmx_data.rng_size[j] = 0;
254     gmmxring.nring_atoms[i] = 0;
255     for (i=0; i < 30; i++)
256     {
257     gmmx_data.ring_atoms[i][j] = 0 ;
258     tmp.tmp_ratoms[i][j] = 0 ;
259     }
260     gmmx_data.clo_distance[j] = 0;
261     tmp.tmp_clo_distance[j] = 0 ;
262     gmmx_data.clo_ang[0][j] = 0;
263     gmmx_data.clo_ang[1][j] = 0;
264     tmp.tmp_clo_ang[0][j] = 0;
265     tmp.tmp_clo_ang[1][j] = 0;
266     gmmx_data.ring_resolution[j] = 0.0;
267     tmp.tmp_ring_resolution[j] = 0.0;
268     gmmx_data.clo_bond[0][j] = 0 ;
269     gmmx_data.clo_bond[1][j] = 0 ;
270     gmmx_data.clo_bond[2][j] = 0 ;
271     gmmx_data.clo_bond[3][j] = 0 ;
272     gmmx_data.clo_bond[4][j] = 0 ;
273     gmmx_data.clo_bond[5][j] = 0 ;
274     tmp.tmp_clo_bond[0][j] = 0 ;
275     tmp.tmp_clo_bond[1][j] = 0 ;
276     tmp.tmp_clo_bond[2][j] = 0 ;
277     tmp.tmp_clo_bond[3][j] = 0 ;
278     }
279     }
280     stmp = rand();
281     if(stmp >= 10000)
282     stmp = stmp/10;
283     if(stmp >= 10000)
284     stmp = stmp/10;
285     if(stmp >= 10000)
286     stmp = stmp/10;
287     if(stmp >= 10000)
288     stmp = stmp/10;
289     if(stmp >= 10000)
290     stmp = stmp/10;
291     if(stmp >= 10000)
292     stmp = stmp/10;
293     if(stmp >= 10000)
294     stmp = stmp/10;
295     if(stmp >= 10000)
296     stmp = stmp/10;
297     sprintf(gmmx_data.jobname,"gmx%d",stmp);
298     strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
299    
300     gmmx_data.method = 3;
301     gmmx_data.iseed = 71277;
302     gmmx_data.its = 990;
303     gmmx_data.lnant = TRUE;
304     gmmx_data.hybrid = TRUE;
305     gmmx_data.ecut = TRUE;
306     gmmx_data.chig = TRUE;
307     gmmx_data.bad15 = TRUE;
308     gmmx_data.nopi1 = FALSE;
309     for (i=1; i <= natom; i++)
310     {
311     if (atom[i].flags & mask)
312     {
313     gmmx_data.nopi1 = TRUE;
314     break;
315     }
316     }
317     gmmx_data.hbond = TRUE;
318     gmmx_data.heat = FALSE;
319    
320     gmmx_data.comp = 0;
321     gmmx_data.qpmr = FALSE;
322     gmmx_data.qdist = FALSE;
323     gmmx_data.qang = FALSE;
324     gmmx_data.qdihed = FALSE;
325     gmmx_data.npmr = 0;
326     gmmx_data.ndist = 0;
327     gmmx_data.nang = 0;
328     gmmx_data.ndihed = 0;
329    
330     gmmx_data.nrings = 0;
331     gmmx_data.nbonds = 0;
332     gmmx_data.ewindow = 3.5;
333     gmmx_data.ewindow2 = 3.0;
334     gmmx_data.boltz = 300.0;
335     gmmx_data.ecutoff = 0.100;
336     gmmx_data.bad15_cutoff = 3.00;
337    
338     gmmx_data.nsrms = 0;
339     gmmx_data.nsrmsa = 0;
340     gmmx_data.comp_method = 0;
341     gmmx_data.comp_method2 = 1;
342     gmmx_data.ermsa = 0.100;
343     gmmx_data.crmsa = 0.250;
344     gmmx_data.restart = FALSE;
345     gmmx_data.include_file = FALSE;
346    
347     nheav = 0;
348     nhyd = 0;
349     for (i=1; i <= natom; i++)
350     {
351     if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
352     nheav++;
353     else
354     nhyd++;
355     }
356    
357     gmmx_data.kstop = 5;
358     gmmx_data.kmin = 5*nheav;
359     gmmx_data.kdup = 5*nheav/2;
360     if (gmmx_data.kdup > 50)
361     gmmx_data.kdup = 50;
362     gmmx_data.max_search = 100000;
363    
364     gmmx_data.nrings = 0;
365     for (i=0; i < 4; i++)
366     gmmx_data.nring_bonds[i] = 0;
367     for (i=0; i < 30; i++)
368     {
369     gmmxring.nring_atoms[i] = 0;
370     for (j=0; j < 30; j++)
371     gmmxring.ring_atoms[i][j]=0;
372     for (j=0; j < 4; j++)
373     for (k=0; k < 3; k++)
374     gmmx_data.ring_bond_data[j][i][k] = 0;
375     }
376    
377     }
378     // =======================================
379     int close_contact(int mode)
380     {
381     int i, j;
382     double dx,dy,dz;
383    
384     for (i=1; i < natom; i++)
385     {
386     for (j=i+1; j <= natom; j++)
387     {
388     dx = fabs(atom[i].x - atom[j].x);
389     dy = fabs(atom[i].y - atom[j].y);
390     dz = fabs(atom[i].z - atom[j].z);
391     if ( (dx + dy + dz) < 0.1)
392     {
393     return FALSE;
394     }
395     }
396     }
397     return TRUE;
398     }