ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 20
Committed: Mon Jun 30 22:00:43 2008 UTC (11 years, 4 months ago) by tjod
Original Path: trunk/smi23d/src/mengine/gmmx_run.c
File size: 14385 byte(s)
Log Message:
Add -a switch to cause H atoms to be added, deleting existing ones.
With -a, mengine behaves as the original.
Formal charges are not always preserved, by the original or
the modified mengine with or without the -a switch.
It appears that negative charges are preserved (benzoic acid).
A positive charge on a N atom (an amine) is not preserved.
Charges are recognized in CHG records and in atom records.
Charges are output only in atom records.

Line File contents
1 #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*,FILE*,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 *paramfile,
165 FILE *infile,
166 int flags)
167 {
168 char line[121];
169 int i,ncount,nret,j,k,jj,nc;
170 int icount;
171 int ngood, nbparam, nbatom,nbcontact;
172 int bcontact = FALSE;
173 int batom = FALSE;
174 float logp;
175 //char *std_file;
176 char bad_atom[80],bad_contact[80],bad_param[80];
177 clock_t start_time, end_time;
178 FILE *logfile;
179
180 // check filetype
181 ncount = 0;
182 input_type = 0;
183
184 /*
185 std_file = strdup(out_file);
186 strcpy(Savebox.fname, out_file);
187 */
188
189 // build filenames
190 strcpy(bad_atom,"_batom.sdf");
191 strcpy(bad_contact,"_bcontact.sdf");
192 strcpy(bad_param,"_bparam.sdf");
193
194 // open logfile
195 logfile = stderr;
196
197 files.nfiles = 0;
198 /*
199 // assume SDF input
200 while (FetchRecord(infile, line) )
201 {
202 if (strcmp(line,"$$$$") == 0)
203 files.nfiles++;
204 }
205 fclose(infile);
206 */
207 //
208 minim_control.field = MMFF94;
209 strcpy(line,paramfile);
210 zero_data();
211 read_datafiles(line);
212 gmmx.run = TRUE;
213 minim_control.method = 1; // just do first deriv minimization
214 icount = 0;
215 //
216 ngood = nbparam = nbatom = nbcontact = 0;
217 while (1)
218 // for (i=1; i <= files.nfiles; i++)
219 {
220 initialize();
221 nret = rd_sdf(infile);
222 if (-1 == nret) return ;
223
224 if (natom == 1)
225 {
226 fprintf(logfile,"Strnum: %d has only one atom\n",i);
227 goto L_10;
228 }
229 // if (i%20 == 0)
230 // printf("Strnum: %d\n",files.nfiles);
231 if (nret == TRUE)
232 {
233 nc = close_contact(1);
234 if (nc == FALSE)
235 {
236 fprintf(logfile,"Strnum: %d CID: %s fails on bad contact\n",files.nfiles,Struct_Title);
237 strcpy(Savebox.fname,bad_contact);
238 sprintf(Struct_Title,"Bad Contact %d\n",i);
239 files.append = TRUE;
240 if (bcontact == FALSE)
241 {
242 files.append = FALSE;
243 bcontact = TRUE;
244 }
245 nbcontact++;
246 write_sdf(flags);
247 //strcpy(Savebox.fname,std_file);
248 goto L_10;
249 }
250 ncount = natom;
251 for (j=1; j<= ncount; j++)
252 {
253 jj = 0;
254 if (atom[j].atomnum == 1) // bridging hydrogens
255 {
256 if (atom[j].iat[0] != 0 && atom[j].iat[1] != 0)
257 goto L_10;
258 }
259 }
260 // TJO
261 if (flags & DO_ADDH) {
262 hdel(0);
263 hadd();
264 }
265 type();
266 // TJO
267 if (flags & DO_ADDH) {
268 hdel(0);
269 }
270 start_time = clock();
271 nret = setup_calculation();
272 if (nret == TRUE)
273 minimize();
274 end_calculation();
275
276 // minim_control.method = 3;
277 // TJO
278 if (flags & DO_ADDH) {
279 hdel(0);
280 hadd();
281 }
282 type_mmx();
283
284 nret = setup_calculation();
285 if (nret == TRUE)
286 {
287 minimize();
288 files.append = TRUE;
289 if ( icount == 0)
290 {
291 files.append = FALSE;
292 icount++;
293 }
294 // compute xlogp
295 if (energies.total < 1000.0)
296 {
297 if (flags & DO_XLOGP) {
298 if (VERBOSE) fprintf(stderr, " Evaluating XlogP\n");
299 xlogp(&logp);
300 logp_calc.logp = logp;
301 }
302 // compute dipole moment
303 dipolemom.xdipole = 0.0;
304 dipolemom.ydipole = 0.0;
305 dipolemom.zdipole = 0.0;
306 if (flags & DO_DIPOLE) {
307 if (VERBOSE) fprintf(stderr, " Evaluating dipole moment\n");
308 charge_dipole();
309 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
310 dipolemom.ydipole*dipolemom.ydipole +
311 dipolemom.zdipole*dipolemom.zdipole);
312 }
313 // compute vib
314 if (flags & DO_VIBRATION) {
315 if (VERBOSE) fprintf(stderr, " Evaluating vibrational parameters\n");
316 vibrate();
317 }
318 write_sdf(flags);
319 ++ngood;
320 ++files.nfiles;
321 }
322 } else
323 {
324 strcpy(Savebox.fname,bad_param);
325 files.append = TRUE;
326 write_sdf(flags);
327 nbparam++;
328 //strcpy(Savebox.fname,std_file);
329 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
330 if (VERBOSE)
331 fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
332 }
333 end_calculation();
334 end_time = clock();
335 fprintf(logfile,"Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
336 natom,(end_time-start_time)/CLOCKS_PER_SEC);
337 if (VERBOSE) {
338 fprintf(stdout,
339 "Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
340 natom,(end_time-start_time)/CLOCKS_PER_SEC);
341 }
342 } else // read_sdf failed
343 {
344 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
345 strcpy(Savebox.fname,bad_atom);
346 files.append = TRUE;
347 if (batom == FALSE)
348 {
349 batom = TRUE;
350 files.append = FALSE;
351 }
352 write_sdf(flags);
353 nbatom++;
354 //strcpy(Savebox.fname,std_file);
355 }
356 L_10:
357 continue;
358 }
359 fclose(infile);
360 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
361 if (VERBOSE) {
362 fprintf(stdout,
363 "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
364 ngood,nbparam,nbatom,nbcontact);
365 }
366
367 fclose(logfile);
368 }
369 // ================================================
370 // ==============================
371 void initialize_gmmx()
372 {
373 int i, j, k, nheav, nhyd,stmp;
374 long int mask;
375
376 mask = (1L << 0);
377 tmp.tmp_nrings = 0;
378 gmmxring.nrings = 0;
379 for (i=0; i < 4; i++)
380 {
381 tmp.tmp_natoms[i] = 0;
382 for (j=0; j < 4; j++)
383 {
384 gmmx_data.rng_size[j] = 0;
385 gmmxring.nring_atoms[i] = 0;
386 for (i=0; i < 30; i++)
387 {
388 gmmx_data.ring_atoms[i][j] = 0 ;
389 tmp.tmp_ratoms[i][j] = 0 ;
390 }
391 gmmx_data.clo_distance[j] = 0;
392 tmp.tmp_clo_distance[j] = 0 ;
393 gmmx_data.clo_ang[0][j] = 0;
394 gmmx_data.clo_ang[1][j] = 0;
395 tmp.tmp_clo_ang[0][j] = 0;
396 tmp.tmp_clo_ang[1][j] = 0;
397 gmmx_data.ring_resolution[j] = 0.0;
398 tmp.tmp_ring_resolution[j] = 0.0;
399 gmmx_data.clo_bond[0][j] = 0 ;
400 gmmx_data.clo_bond[1][j] = 0 ;
401 gmmx_data.clo_bond[2][j] = 0 ;
402 gmmx_data.clo_bond[3][j] = 0 ;
403 gmmx_data.clo_bond[4][j] = 0 ;
404 gmmx_data.clo_bond[5][j] = 0 ;
405 tmp.tmp_clo_bond[0][j] = 0 ;
406 tmp.tmp_clo_bond[1][j] = 0 ;
407 tmp.tmp_clo_bond[2][j] = 0 ;
408 tmp.tmp_clo_bond[3][j] = 0 ;
409 }
410 }
411 stmp = rand();
412 if(stmp >= 10000)
413 stmp = stmp/10;
414 if(stmp >= 10000)
415 stmp = stmp/10;
416 if(stmp >= 10000)
417 stmp = stmp/10;
418 if(stmp >= 10000)
419 stmp = stmp/10;
420 if(stmp >= 10000)
421 stmp = stmp/10;
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 sprintf(gmmx_data.jobname,"gmx%d",stmp);
429 strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
430
431 gmmx_data.method = 3;
432 gmmx_data.iseed = 71277;
433 gmmx_data.its = 990;
434 gmmx_data.lnant = TRUE;
435 gmmx_data.hybrid = TRUE;
436 gmmx_data.ecut = TRUE;
437 gmmx_data.chig = TRUE;
438 gmmx_data.bad15 = TRUE;
439 gmmx_data.nopi1 = FALSE;
440 for (i=1; i <= natom; i++)
441 {
442 if (atom[i].flags & mask)
443 {
444 gmmx_data.nopi1 = TRUE;
445 break;
446 }
447 }
448 gmmx_data.hbond = TRUE;
449 gmmx_data.heat = FALSE;
450
451 gmmx_data.comp = 0;
452 gmmx_data.qpmr = FALSE;
453 gmmx_data.qdist = FALSE;
454 gmmx_data.qang = FALSE;
455 gmmx_data.qdihed = FALSE;
456 gmmx_data.npmr = 0;
457 gmmx_data.ndist = 0;
458 gmmx_data.nang = 0;
459 gmmx_data.ndihed = 0;
460
461 gmmx_data.nrings = 0;
462 gmmx_data.nbonds = 0;
463 gmmx_data.ewindow = 3.5;
464 gmmx_data.ewindow2 = 3.0;
465 gmmx_data.boltz = 300.0;
466 gmmx_data.ecutoff = 0.100;
467 gmmx_data.bad15_cutoff = 3.00;
468
469 gmmx_data.nsrms = 0;
470 gmmx_data.nsrmsa = 0;
471 gmmx_data.comp_method = 0;
472 gmmx_data.comp_method2 = 1;
473 gmmx_data.ermsa = 0.100;
474 gmmx_data.crmsa = 0.250;
475 gmmx_data.restart = FALSE;
476 gmmx_data.include_file = FALSE;
477
478 nheav = 0;
479 nhyd = 0;
480 for (i=1; i <= natom; i++)
481 {
482 if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
483 nheav++;
484 else
485 nhyd++;
486 }
487
488 gmmx_data.kstop = 5;
489 gmmx_data.kmin = 5*nheav;
490 gmmx_data.kdup = 5*nheav/2;
491 if (gmmx_data.kdup > 50)
492 gmmx_data.kdup = 50;
493 gmmx_data.max_search = 100000;
494
495 gmmx_data.nrings = 0;
496 for (i=0; i < 4; i++)
497 gmmx_data.nring_bonds[i] = 0;
498 for (i=0; i < 30; i++)
499 {
500 gmmxring.nring_atoms[i] = 0;
501 for (j=0; j < 30; j++)
502 gmmxring.ring_atoms[i][j]=0;
503 for (j=0; j < 4; j++)
504 for (k=0; k < 3; k++)
505 gmmx_data.ring_bond_data[j][i][k] = 0;
506 }
507
508 }
509 /* =============================================== */
510 // get hybridization of atom - tetrahedral = 1, planar = 2, linear = 3
511 //
512 int get_hybrid(int ia)
513 {
514 int itype;
515
516 itype = atom[ia].mmx_type;
517 if (itype == 2 || itype == 3 || itype == 7 || itype == 9 || itype == 25 ||
518 itype == 26 || itype == 29 || itype == 30 || itype == 37 || itype == 38 ||
519 itype == 40 || itype == 57 )
520 return 2;
521 else if (itype == 4 || itype == 10)
522 return 3;
523 else
524 return 1;
525 }
526 // ==================================
527 int check_ring1(int ia)
528 {
529 if (is_ring61(ia)) return TRUE;
530 if (is_ring51(ia)) return TRUE;
531 if (is_ring31(ia)) return TRUE;
532 if (is_ring41(ia)) return TRUE;
533 return FALSE;
534 }
535 // =======================================
536 int close_contact(int mode)
537 {
538 int i, j;
539 double dx,dy,dz;
540
541 for (i=1; i < natom; i++)
542 {
543 for (j=i+1; j <= natom; j++)
544 {
545 dx = fabs(atom[i].x - atom[j].x);
546 dy = fabs(atom[i].y - atom[j].y);
547 dz = fabs(atom[i].z - atom[j].z);
548 if ( (dx + dy + dz) < 0.1)
549 {
550 return FALSE;
551 }
552 }
553 }
554 return TRUE;
555 }