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