ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 75
Committed: Mon Dec 15 18:03:29 2008 UTC (13 years, 5 months ago) by gilbertke
File size: 9839 byte(s)
Log Message:
fixes in typeing to pass MMFF validation suite
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 75 // minim_values.iprint = TRUE;
156 tjod 3 nret = setup_calculation();
157 gilbertke 75 // for(i=1; i<= natom; i++)
158     //printf("Atom: %d %d %x\n",i,atom[i].type,atom[i].flags);
159     //if (minim_values.iprint == TRUE) exit(0);
160 tjod 3 if (nret == TRUE)
161     {
162     minimize();
163     // compute xlogp
164     if (energies.total < 1000.0)
165     {
166     if (flags & DO_XLOGP) {
167 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating XlogP\n");
168 tjod 3 xlogp(&logp);
169     logp_calc.logp = logp;
170     }
171     // compute dipole moment
172     dipolemom.xdipole = 0.0;
173     dipolemom.ydipole = 0.0;
174     dipolemom.zdipole = 0.0;
175     if (flags & DO_DIPOLE) {
176 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating dipole moment\n");
177 tjod 3 charge_dipole();
178     dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
179     dipolemom.ydipole*dipolemom.ydipole +
180     dipolemom.zdipole*dipolemom.zdipole);
181     }
182     // compute vib
183     if (flags & DO_VIBRATION) {
184 wdelano 58 if (VERBOSE) fprintf(pcmlogfile, " Evaluating vibrational parameters\n");
185 tjod 3 vibrate();
186     }
187     write_sdf(flags);
188 tjod 15 ++ngood;
189     ++files.nfiles;
190 tjod 3 }
191     } else
192     {
193     strcpy(Savebox.fname,bad_param);
194     files.append = TRUE;
195 tjod 33 //write_sdf(flags);
196 tjod 3 nbparam++;
197 tjod 15 //strcpy(Savebox.fname,std_file);
198     fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
199 tjod 3 if (VERBOSE)
200     fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
201     }
202     end_calculation();
203     if (VERBOSE) {
204     fprintf(stdout,
205 gilbertke 63 "Strnum: %d %s natom %d \n",files.nfiles,Struct_Title, natom);
206 tjod 3 }
207     } else // read_sdf failed
208     {
209 tjod 15 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
210 tjod 3 strcpy(Savebox.fname,bad_atom);
211     files.append = TRUE;
212     if (batom == FALSE)
213     {
214     batom = TRUE;
215     files.append = FALSE;
216     }
217 tjod 33 //write_sdf(flags);
218 tjod 3 nbatom++;
219 tjod 15 //strcpy(Savebox.fname,std_file);
220 tjod 3 }
221     L_10:
222     continue;
223     }
224 wdelano 58 // fclose(infile);
225 tjod 3 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
226     if (VERBOSE) {
227     fprintf(stdout,
228     "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
229     ngood,nbparam,nbatom,nbcontact);
230     }
231    
232 wdelano 58 // fclose(logfile);
233 tjod 3 }
234     // ================================================
235     // ==============================
236     void initialize_gmmx()
237     {
238     int i, j, k, nheav, nhyd,stmp;
239     long int mask;
240    
241     mask = (1L << 0);
242     tmp.tmp_nrings = 0;
243     gmmxring.nrings = 0;
244     for (i=0; i < 4; i++)
245     {
246     tmp.tmp_natoms[i] = 0;
247     for (j=0; j < 4; j++)
248     {
249     gmmx_data.rng_size[j] = 0;
250     gmmxring.nring_atoms[i] = 0;
251     for (i=0; i < 30; i++)
252     {
253     gmmx_data.ring_atoms[i][j] = 0 ;
254     tmp.tmp_ratoms[i][j] = 0 ;
255     }
256     gmmx_data.clo_distance[j] = 0;
257     tmp.tmp_clo_distance[j] = 0 ;
258     gmmx_data.clo_ang[0][j] = 0;
259     gmmx_data.clo_ang[1][j] = 0;
260     tmp.tmp_clo_ang[0][j] = 0;
261     tmp.tmp_clo_ang[1][j] = 0;
262     gmmx_data.ring_resolution[j] = 0.0;
263     tmp.tmp_ring_resolution[j] = 0.0;
264     gmmx_data.clo_bond[0][j] = 0 ;
265     gmmx_data.clo_bond[1][j] = 0 ;
266     gmmx_data.clo_bond[2][j] = 0 ;
267     gmmx_data.clo_bond[3][j] = 0 ;
268     gmmx_data.clo_bond[4][j] = 0 ;
269     gmmx_data.clo_bond[5][j] = 0 ;
270     tmp.tmp_clo_bond[0][j] = 0 ;
271     tmp.tmp_clo_bond[1][j] = 0 ;
272     tmp.tmp_clo_bond[2][j] = 0 ;
273     tmp.tmp_clo_bond[3][j] = 0 ;
274     }
275     }
276     stmp = rand();
277     if(stmp >= 10000)
278     stmp = stmp/10;
279     if(stmp >= 10000)
280     stmp = stmp/10;
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     sprintf(gmmx_data.jobname,"gmx%d",stmp);
294     strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
295    
296     gmmx_data.method = 3;
297     gmmx_data.iseed = 71277;
298     gmmx_data.its = 990;
299     gmmx_data.lnant = TRUE;
300     gmmx_data.hybrid = TRUE;
301     gmmx_data.ecut = TRUE;
302     gmmx_data.chig = TRUE;
303     gmmx_data.bad15 = TRUE;
304     gmmx_data.nopi1 = FALSE;
305     for (i=1; i <= natom; i++)
306     {
307     if (atom[i].flags & mask)
308     {
309     gmmx_data.nopi1 = TRUE;
310     break;
311     }
312     }
313     gmmx_data.hbond = TRUE;
314     gmmx_data.heat = FALSE;
315    
316     gmmx_data.comp = 0;
317     gmmx_data.qpmr = FALSE;
318     gmmx_data.qdist = FALSE;
319     gmmx_data.qang = FALSE;
320     gmmx_data.qdihed = FALSE;
321     gmmx_data.npmr = 0;
322     gmmx_data.ndist = 0;
323     gmmx_data.nang = 0;
324     gmmx_data.ndihed = 0;
325    
326     gmmx_data.nrings = 0;
327     gmmx_data.nbonds = 0;
328     gmmx_data.ewindow = 3.5;
329     gmmx_data.ewindow2 = 3.0;
330     gmmx_data.boltz = 300.0;
331     gmmx_data.ecutoff = 0.100;
332     gmmx_data.bad15_cutoff = 3.00;
333    
334     gmmx_data.nsrms = 0;
335     gmmx_data.nsrmsa = 0;
336     gmmx_data.comp_method = 0;
337     gmmx_data.comp_method2 = 1;
338     gmmx_data.ermsa = 0.100;
339     gmmx_data.crmsa = 0.250;
340     gmmx_data.restart = FALSE;
341     gmmx_data.include_file = FALSE;
342    
343     nheav = 0;
344     nhyd = 0;
345     for (i=1; i <= natom; i++)
346     {
347     if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
348     nheav++;
349     else
350     nhyd++;
351     }
352    
353     gmmx_data.kstop = 5;
354     gmmx_data.kmin = 5*nheav;
355     gmmx_data.kdup = 5*nheav/2;
356     if (gmmx_data.kdup > 50)
357     gmmx_data.kdup = 50;
358     gmmx_data.max_search = 100000;
359    
360     gmmx_data.nrings = 0;
361     for (i=0; i < 4; i++)
362     gmmx_data.nring_bonds[i] = 0;
363     for (i=0; i < 30; i++)
364     {
365     gmmxring.nring_atoms[i] = 0;
366     for (j=0; j < 30; j++)
367     gmmxring.ring_atoms[i][j]=0;
368     for (j=0; j < 4; j++)
369     for (k=0; k < 3; k++)
370     gmmx_data.ring_bond_data[j][i][k] = 0;
371     }
372    
373     }
374     // =======================================
375     int close_contact(int mode)
376     {
377     int i, j;
378     double dx,dy,dz;
379    
380     for (i=1; i < natom; i++)
381     {
382     for (j=i+1; j <= natom; j++)
383     {
384     dx = fabs(atom[i].x - atom[j].x);
385     dy = fabs(atom[i].y - atom[j].y);
386     dz = fabs(atom[i].z - atom[j].z);
387     if ( (dx + dy + dz) < 0.1)
388     {
389     return FALSE;
390     }
391     }
392     }
393     return TRUE;
394     }