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 (10 years, 11 months ago) by gilbertke
File size: 9849 byte(s)
Log Message:
added aromatic bonds in read_sdf removed TJO bailout
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 #include "dipmom.h"
13 #include "vibrate.h"
14
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
53 int setup_calculation(void);
54 void end_calculation(void);
55 void minimize(void);
56 static int close_contact(int);
57
58 void initialize(void);
59
60 void type(void);
61
62 void read_gmmxinp(char*,FILE*,int);
63 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 void read_gmmxinp(char *paramfile,
74 FILE *infile,
75 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 //char *std_file;
85 char bad_atom[80],bad_contact[80],bad_param[80];
86 FILE *logfile;
87
88 // check filetype
89 ncount = 0;
90 input_type = 0;
91
92 /*
93 std_file = strdup(out_file);
94 strcpy(Savebox.fname, out_file);
95 */
96
97 // build filenames
98 strcpy(bad_atom,"_batom.sdf");
99 strcpy(bad_contact,"_bcontact.sdf");
100 strcpy(bad_param,"_bparam.sdf");
101
102 // open logfile
103 logfile = pcmlogfile;
104
105 files.nfiles = 0;
106 /*
107 // assume SDF input
108 while (FetchRecord(infile, line) )
109 {
110 if (strcmp(line,"$$$$") == 0)
111 files.nfiles++;
112 }
113 fclose(infile);
114 */
115 //
116 minim_control.field = MMFF94;
117 //strcpy(line,paramfile);
118 zero_data();
119 read_datafiles(paramfile);
120 gmmx.run = TRUE;
121 minim_control.method = 1; // just do first deriv minimization
122 icount = 0;
123 //
124 ngood = nbparam = nbatom = nbcontact = 0;
125 i = 0;
126 while (1)
127 // for (i=1; i <= files.nfiles; i++)
128 {
129 ++i;
130 initialize();
131 nret = rd_sdf(infile);
132 if (-1 == nret) return ;
133 if (FALSE == nret)
134 {
135 fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
136 goto L_10;
137 natom = 0;
138 }
139
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 // TJO
148 if (flags & DO_ADDH)
149 {
150 hdel(0);
151 hadd();
152 }
153 type();
154
155 // minim_values.iprint = TRUE;
156 nret = setup_calculation();
157 // if (minim_values.iprint == TRUE) exit(0);
158 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 if (VERBOSE) fprintf(pcmlogfile, " Evaluating XlogP\n");
172 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 if (VERBOSE) fprintf(pcmlogfile, " Evaluating dipole moment\n");
181 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 if (VERBOSE) fprintf(pcmlogfile, " Evaluating vibrational parameters\n");
189 vibrate();
190 }
191 write_sdf(flags);
192 ++ngood;
193 ++files.nfiles;
194 }
195 } else
196 {
197 strcpy(Savebox.fname,bad_param);
198 files.append = TRUE;
199 //write_sdf(flags);
200 nbparam++;
201 //strcpy(Savebox.fname,std_file);
202 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
203 if (VERBOSE)
204 fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
205 }
206 end_calculation();
207 if (VERBOSE) {
208 fprintf(stdout,
209 "Strnum: %d %s natom %d \n",files.nfiles,Struct_Title, natom);
210 }
211 } else // read_sdf failed
212 {
213 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
214 strcpy(Savebox.fname,bad_atom);
215 files.append = TRUE;
216 if (batom == FALSE)
217 {
218 batom = TRUE;
219 files.append = FALSE;
220 }
221 //write_sdf(flags);
222 nbatom++;
223 //strcpy(Savebox.fname,std_file);
224 }
225 L_10:
226 continue;
227 }
228 // fclose(infile);
229 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 // fclose(logfile);
237 }
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 }