ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 48
Committed: Thu Jul 31 16:34:52 2008 UTC (12 years, 10 months ago) by tjod
File size: 14648 byte(s)
Log Message:
Add modification notice to files modified for Freemol
per mengine/smi23d license requirements

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    
13     #include <time.h>
14     #include <errno.h>
15    
16     struct t_logp {
17     float logp;
18     } logp_calc;
19     struct t_vibdata {
20     char ptgrp[4];
21     float mom_ix,mom_iy,mom_iz;
22     float etot,htot,stot,gtot,cptot;
23     } vibdata;
24    
25     EXTERN struct t_optimize {
26     int param_avail, converge;
27     float initial_energy, final_energy, initial_heat, final_heat;
28     } optimize_data;
29    
30     EXTERN struct t_files {
31     int nfiles, append, batch, icurrent, ibatno;
32     } files;
33     EXTERN struct t_dipolemom {
34     double total, xdipole, ydipole, zdipole;
35     } dipolemom;
36    
37     EXTERN struct t_hfof {
38     float hfo, si;
39     int iunk;
40     } hfof;
41     EXTERN struct t_pcmfile {
42     char string[200];
43     int head;
44     char token[20];
45     int state;
46     unsigned int nocaps;
47     } pcmfile;
48     EXTERN struct t_minim_values {
49     int iprint, ndc, nconst;
50     float dielc;
51     } minim_values;
52     EXTERN struct t_minim_control {
53     int type, method, field, added_const;
54     char added_path[256],added_name[256];
55     } minim_control;
56     EXTERN struct t_rotbond {
57     int nbonds, bond[200][2], bond_res[200];
58     int incl_amide, incl_alkene;
59     } rotbond;
60    
61     struct t_tmp {
62     int tmp_nrings, tmp_natoms[4], tmp_ratoms[30][4], tmp_clo_bond[4][4];
63     float tmp_clo_ang[2][4], tmp_clo_distance[4], tmp_ring_resolution[4];
64     } tmp;
65    
66     int input_type;
67    
68     int nerr;
69     int gmmx_abort;
70    
71     int rd_sdf(FILE *);
72     int check_ring1(int);
73     int is_ring31(int);
74     int is_ring41(int);
75     int is_ring51(int);
76     int is_ring61(int);
77     int get_hybrid(int);
78     void set_atomtype(int,int,int,int,int,int);
79     void sort_file(char *, char *);
80     void gettoken(void);
81     int is_bond(int,int);
82     void read_atomtypes(int);
83     //double compute_rms(int, int **, double [150][3], double [150][3]);
84     //void fndumt( int, int **, double [150][3], double [150][3] );
85     //int fast_compare(double [150][3], double [150][3]);
86     //void slow_compare(int, int **, double [150][3], double [150][3]);
87     double compute_rms(int, int **, double **, double **);
88     void fndumt( int, int **, double **, double ** );
89     int fast_compare(double **, double **);
90     void slow_compare(int, int **, double **, double **);
91     void message_alert(char *, char *);
92     void gradient(void);
93     void eigen(int norder,double (*)[3],double *, double (*)[3]);
94     int check_structure(void);
95     int check_id(void);
96     double distance(int,int);
97     void rotb(int,int,int,int, float);
98     void make_bond(int , int , int );
99     void deletebond(int, int);
100     int setup_calculation(void);
101     void end_calculation(void);
102     void hster(int, int *);
103     int lepimer(int);
104     double dihdrl(int,int,int,int);
105     void mvatom(void);
106     void mvbond(void);
107     void minimize(void);
108     int findcd(int);
109     int findtor(int);
110     void InitialTransform(void);
111     int bad15(int);
112     int close_contact(int);
113     void gmmx2(void);
114     int check_stop(void);
115     void hbondreset(void);
116     void setup_output(char *);
117     void pimarkselection(void);
118     void mark_hbond(void);
119     void pireset(void);
120     void pcmfout(int);
121     void center(int);
122     void pcmfin(int,int);
123     void initialize(void);
124     void generate_bonds(void);
125     int bstereo(int);
126     void type(void);
127     void mmxsub(int);
128     void set_atom_color(void);
129     void eheat(void);
130     void charge_dipole(void);
131     void mmp22mod(int,int);
132     void check_numfile(int);
133     void dipole_dipole(void);
134     void set_active(void);
135     double energy(void);
136     void write_restart(void);
137     int read_restart(void);
138     void rdfile(int,int);
139     void montc(void);
140     void include_min(void);
141     void jacobi(int,int,double **,double *,double **,double *,double *);
142     void quatfit(int,int **,double **,double **);
143     FILE * fopen_path ( char * , char * , char * ) ;
144 tjod 16 void read_gmmxinp(char*,FILE*,int);
145 tjod 3 void run_gmmx(void);
146     int rmsck(int, int *);
147     void write_gmmx(void);
148     int read_sdf(int,int);
149     float ran1(int *);
150     void hdel(int);
151     void hadd(void);
152     void hcoord(void);
153     void inesc(char *);
154     void find_bonds(void);
155     int FetchRecord(FILE *, char *);
156     void deleteatom(int);
157     void write_sdf(int);
158     void zero_data(void);
159     void type_mmx(void);
160     void read_datafiles(char *);
161     void xlogp(float *);
162     void vibrate(void);
163    
164     // =====================================
165 tjod 16 void read_gmmxinp(char *paramfile,
166 tjod 15 FILE *infile,
167 tjod 3 int flags)
168     {
169     char line[121];
170     int i,ncount,nret,j,k,jj,nc;
171     int icount;
172     int ngood, nbparam, nbatom,nbcontact;
173     int bcontact = FALSE;
174     int batom = FALSE;
175     float logp;
176 tjod 15 //char *std_file;
177 tjod 3 char bad_atom[80],bad_contact[80],bad_param[80];
178     clock_t start_time, end_time;
179 tjod 15 FILE *logfile;
180 tjod 3
181     // check filetype
182     ncount = 0;
183     input_type = 0;
184    
185 tjod 15 /*
186 tjod 3 std_file = strdup(out_file);
187     strcpy(Savebox.fname, out_file);
188 tjod 15 */
189 tjod 3
190     // build filenames
191 tjod 15 strcpy(bad_atom,"_batom.sdf");
192     strcpy(bad_contact,"_bcontact.sdf");
193     strcpy(bad_param,"_bparam.sdf");
194 tjod 3
195     // open logfile
196 tjod 16 logfile = stderr;
197 tjod 3
198     files.nfiles = 0;
199 tjod 15 /*
200 tjod 3 // assume SDF input
201     while (FetchRecord(infile, line) )
202     {
203     if (strcmp(line,"$$$$") == 0)
204     files.nfiles++;
205     }
206     fclose(infile);
207 tjod 15 */
208 tjod 3 //
209     minim_control.field = MMFF94;
210 tjod 25 //strcpy(line,paramfile);
211 tjod 3 zero_data();
212 tjod 25 read_datafiles(paramfile);
213 tjod 3 gmmx.run = TRUE;
214     minim_control.method = 1; // just do first deriv minimization
215     icount = 0;
216     //
217     ngood = nbparam = nbatom = nbcontact = 0;
218 tjod 31 i = 0;
219 tjod 15 while (1)
220     // for (i=1; i <= files.nfiles; i++)
221 tjod 3 {
222 tjod 31 ++i;
223 tjod 3 initialize();
224     nret = rd_sdf(infile);
225 tjod 15 if (-1 == nret) return ;
226 tjod 31 if (FALSE == nret) {
227     fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
228     goto L_10;
229     natom = 0;
230     }
231 tjod 3
232     if (natom == 1)
233     {
234     fprintf(logfile,"Strnum: %d has only one atom\n",i);
235     goto L_10;
236     }
237     // if (i%20 == 0)
238 tjod 15 // printf("Strnum: %d\n",files.nfiles);
239 tjod 3 if (nret == TRUE)
240     {
241     nc = close_contact(1);
242     if (nc == FALSE)
243     {
244 tjod 15 fprintf(logfile,"Strnum: %d CID: %s fails on bad contact\n",files.nfiles,Struct_Title);
245 tjod 3 strcpy(Savebox.fname,bad_contact);
246     sprintf(Struct_Title,"Bad Contact %d\n",i);
247     files.append = TRUE;
248     if (bcontact == FALSE)
249     {
250     files.append = FALSE;
251     bcontact = TRUE;
252     }
253     nbcontact++;
254 tjod 33 //write_sdf(flags);
255 tjod 15 //strcpy(Savebox.fname,std_file);
256 tjod 3 goto L_10;
257     }
258     ncount = natom;
259     for (j=1; j<= ncount; j++)
260     {
261     jj = 0;
262     if (atom[j].atomnum == 1) // bridging hydrogens
263     {
264     if (atom[j].iat[0] != 0 && atom[j].iat[1] != 0)
265     goto L_10;
266     }
267     }
268 tjod 20 // TJO
269     if (flags & DO_ADDH) {
270     hdel(0);
271     hadd();
272     }
273 tjod 3 type();
274 tjod 20 // TJO
275     if (flags & DO_ADDH) {
276     hdel(0);
277     }
278 tjod 3 start_time = clock();
279     nret = setup_calculation();
280     if (nret == TRUE)
281     minimize();
282     end_calculation();
283    
284     // minim_control.method = 3;
285 tjod 20 // TJO
286     if (flags & DO_ADDH) {
287     hdel(0);
288     hadd();
289     }
290 tjod 3 type_mmx();
291    
292     nret = setup_calculation();
293     if (nret == TRUE)
294     {
295     minimize();
296     files.append = TRUE;
297     if ( icount == 0)
298     {
299     files.append = FALSE;
300     icount++;
301     }
302     // compute xlogp
303     if (energies.total < 1000.0)
304     {
305     if (flags & DO_XLOGP) {
306     if (VERBOSE) fprintf(stderr, " Evaluating XlogP\n");
307     xlogp(&logp);
308     logp_calc.logp = logp;
309     }
310     // compute dipole moment
311     dipolemom.xdipole = 0.0;
312     dipolemom.ydipole = 0.0;
313     dipolemom.zdipole = 0.0;
314     if (flags & DO_DIPOLE) {
315     if (VERBOSE) fprintf(stderr, " Evaluating dipole moment\n");
316     charge_dipole();
317     dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
318     dipolemom.ydipole*dipolemom.ydipole +
319     dipolemom.zdipole*dipolemom.zdipole);
320     }
321     // compute vib
322     if (flags & DO_VIBRATION) {
323     if (VERBOSE) fprintf(stderr, " Evaluating vibrational parameters\n");
324     vibrate();
325     }
326     write_sdf(flags);
327 tjod 15 ++ngood;
328     ++files.nfiles;
329 tjod 3 }
330     } else
331     {
332     strcpy(Savebox.fname,bad_param);
333     files.append = TRUE;
334 tjod 33 //write_sdf(flags);
335 tjod 3 nbparam++;
336 tjod 15 //strcpy(Savebox.fname,std_file);
337     fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,Struct_Title);
338 tjod 3 if (VERBOSE)
339     fprintf(stdout, "Missing parameters - no calc - %s\n",Struct_Title);
340     }
341     end_calculation();
342     end_time = clock();
343 tjod 15 fprintf(logfile,"Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
344 tjod 3 natom,(end_time-start_time)/CLOCKS_PER_SEC);
345     if (VERBOSE) {
346     fprintf(stdout,
347 tjod 15 "Strnum: %d %s natom %d time %lu sec\n",files.nfiles,Struct_Title,
348 tjod 3 natom,(end_time-start_time)/CLOCKS_PER_SEC);
349     }
350     } else // read_sdf failed
351     {
352 tjod 15 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,Struct_Title);
353 tjod 3 strcpy(Savebox.fname,bad_atom);
354     files.append = TRUE;
355     if (batom == FALSE)
356     {
357     batom = TRUE;
358     files.append = FALSE;
359     }
360 tjod 33 //write_sdf(flags);
361 tjod 3 nbatom++;
362 tjod 15 //strcpy(Savebox.fname,std_file);
363 tjod 3 }
364     L_10:
365     continue;
366     }
367     fclose(infile);
368     fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
369     if (VERBOSE) {
370     fprintf(stdout,
371     "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
372     ngood,nbparam,nbatom,nbcontact);
373     }
374    
375     fclose(logfile);
376     }
377     // ================================================
378     // ==============================
379     void initialize_gmmx()
380     {
381     int i, j, k, nheav, nhyd,stmp;
382     long int mask;
383    
384     mask = (1L << 0);
385     tmp.tmp_nrings = 0;
386     gmmxring.nrings = 0;
387     for (i=0; i < 4; i++)
388     {
389     tmp.tmp_natoms[i] = 0;
390     for (j=0; j < 4; j++)
391     {
392     gmmx_data.rng_size[j] = 0;
393     gmmxring.nring_atoms[i] = 0;
394     for (i=0; i < 30; i++)
395     {
396     gmmx_data.ring_atoms[i][j] = 0 ;
397     tmp.tmp_ratoms[i][j] = 0 ;
398     }
399     gmmx_data.clo_distance[j] = 0;
400     tmp.tmp_clo_distance[j] = 0 ;
401     gmmx_data.clo_ang[0][j] = 0;
402     gmmx_data.clo_ang[1][j] = 0;
403     tmp.tmp_clo_ang[0][j] = 0;
404     tmp.tmp_clo_ang[1][j] = 0;
405     gmmx_data.ring_resolution[j] = 0.0;
406     tmp.tmp_ring_resolution[j] = 0.0;
407     gmmx_data.clo_bond[0][j] = 0 ;
408     gmmx_data.clo_bond[1][j] = 0 ;
409     gmmx_data.clo_bond[2][j] = 0 ;
410     gmmx_data.clo_bond[3][j] = 0 ;
411     gmmx_data.clo_bond[4][j] = 0 ;
412     gmmx_data.clo_bond[5][j] = 0 ;
413     tmp.tmp_clo_bond[0][j] = 0 ;
414     tmp.tmp_clo_bond[1][j] = 0 ;
415     tmp.tmp_clo_bond[2][j] = 0 ;
416     tmp.tmp_clo_bond[3][j] = 0 ;
417     }
418     }
419     stmp = rand();
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     if(stmp >= 10000)
429     stmp = stmp/10;
430     if(stmp >= 10000)
431     stmp = stmp/10;
432     if(stmp >= 10000)
433     stmp = stmp/10;
434     if(stmp >= 10000)
435     stmp = stmp/10;
436     sprintf(gmmx_data.jobname,"gmx%d",stmp);
437     strcpy(gmmx_data.gmmx_comment,"GMMX Conf Search");
438    
439     gmmx_data.method = 3;
440     gmmx_data.iseed = 71277;
441     gmmx_data.its = 990;
442     gmmx_data.lnant = TRUE;
443     gmmx_data.hybrid = TRUE;
444     gmmx_data.ecut = TRUE;
445     gmmx_data.chig = TRUE;
446     gmmx_data.bad15 = TRUE;
447     gmmx_data.nopi1 = FALSE;
448     for (i=1; i <= natom; i++)
449     {
450     if (atom[i].flags & mask)
451     {
452     gmmx_data.nopi1 = TRUE;
453     break;
454     }
455     }
456     gmmx_data.hbond = TRUE;
457     gmmx_data.heat = FALSE;
458    
459     gmmx_data.comp = 0;
460     gmmx_data.qpmr = FALSE;
461     gmmx_data.qdist = FALSE;
462     gmmx_data.qang = FALSE;
463     gmmx_data.qdihed = FALSE;
464     gmmx_data.npmr = 0;
465     gmmx_data.ndist = 0;
466     gmmx_data.nang = 0;
467     gmmx_data.ndihed = 0;
468    
469     gmmx_data.nrings = 0;
470     gmmx_data.nbonds = 0;
471     gmmx_data.ewindow = 3.5;
472     gmmx_data.ewindow2 = 3.0;
473     gmmx_data.boltz = 300.0;
474     gmmx_data.ecutoff = 0.100;
475     gmmx_data.bad15_cutoff = 3.00;
476    
477     gmmx_data.nsrms = 0;
478     gmmx_data.nsrmsa = 0;
479     gmmx_data.comp_method = 0;
480     gmmx_data.comp_method2 = 1;
481     gmmx_data.ermsa = 0.100;
482     gmmx_data.crmsa = 0.250;
483     gmmx_data.restart = FALSE;
484     gmmx_data.include_file = FALSE;
485    
486     nheav = 0;
487     nhyd = 0;
488     for (i=1; i <= natom; i++)
489     {
490     if (atom[i].atomnum != 1 && atom[i].atomnum != 2)
491     nheav++;
492     else
493     nhyd++;
494     }
495    
496     gmmx_data.kstop = 5;
497     gmmx_data.kmin = 5*nheav;
498     gmmx_data.kdup = 5*nheav/2;
499     if (gmmx_data.kdup > 50)
500     gmmx_data.kdup = 50;
501     gmmx_data.max_search = 100000;
502    
503     gmmx_data.nrings = 0;
504     for (i=0; i < 4; i++)
505     gmmx_data.nring_bonds[i] = 0;
506     for (i=0; i < 30; i++)
507     {
508     gmmxring.nring_atoms[i] = 0;
509     for (j=0; j < 30; j++)
510     gmmxring.ring_atoms[i][j]=0;
511     for (j=0; j < 4; j++)
512     for (k=0; k < 3; k++)
513     gmmx_data.ring_bond_data[j][i][k] = 0;
514     }
515    
516     }
517     /* =============================================== */
518     // get hybridization of atom - tetrahedral = 1, planar = 2, linear = 3
519     //
520     int get_hybrid(int ia)
521     {
522     int itype;
523    
524     itype = atom[ia].mmx_type;
525     if (itype == 2 || itype == 3 || itype == 7 || itype == 9 || itype == 25 ||
526     itype == 26 || itype == 29 || itype == 30 || itype == 37 || itype == 38 ||
527     itype == 40 || itype == 57 )
528     return 2;
529     else if (itype == 4 || itype == 10)
530     return 3;
531     else
532     return 1;
533     }
534     // ==================================
535     int check_ring1(int ia)
536     {
537     if (is_ring61(ia)) return TRUE;
538     if (is_ring51(ia)) return TRUE;
539     if (is_ring31(ia)) return TRUE;
540     if (is_ring41(ia)) return TRUE;
541     return FALSE;
542     }
543     // =======================================
544     int close_contact(int mode)
545     {
546     int i, j;
547     double dx,dy,dz;
548    
549     for (i=1; i < natom; i++)
550     {
551     for (j=i+1; j <= natom; j++)
552     {
553     dx = fabs(atom[i].x - atom[j].x);
554     dy = fabs(atom[i].y - atom[j].y);
555     dz = fabs(atom[i].z - atom[j].z);
556     if ( (dx + dy + dz) < 0.1)
557     {
558     return FALSE;
559     }
560     }
561     }
562     return TRUE;
563     }