ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/gmmx_run.c
Revision: 127
Committed: Fri Jun 19 18:18:02 2009 UTC (11 years ago) by gilbertke
File size: 6039 byte(s)
Log Message:
ran mmff94.sdf test suite on latest code - fixed errors in rings.c
Line File contents
1 /* NOTICE: this source code file has been modified for use with FreeMOL */
2 #define notKEG_DEBUG
3
4
5 #define EXTERN extern
6 #include "pcwin.h"
7 #include "pcmod.h"
8 #include "utility.h"
9 #include "vibrate.h"
10
11 #include <errno.h>
12
13 struct t_logp {
14 float logp;
15 } logp_calc;
16
17 EXTERN struct t_files {
18 int nfiles, append, batch, icurrent, ibatno;
19 } files;
20
21 EXTERN struct t_minim_values {
22 int iprint, ndc, nconst;
23 float dielc;
24 } minim_values;
25 EXTERN struct t_minim_control {
26 int type, method, field, added_const;
27 char added_path[256],added_name[256];
28 } minim_control;
29
30 int rd_sdf(FILE *);
31 int setup_calculation(void);
32 void end_calculation(void);
33 void minimize(int natom,int *use,double *x,double *y,double *z);
34 static int close_contact(int);
35 void initialize(void);
36 void type(void);
37 void read_gmmxinp(char*,FILE*,int);
38 void hdel(int);
39 void hadd(void);
40 int FetchRecord(FILE *, char *);
41 void write_sdf(int);
42 void zero_data(void);
43 void type_mmx(void);
44 void read_datafiles(char *);
45 float xlogp(int natom,int *atomnum,int **iat,int **bo,long int *flags,int *type);
46 void print_energy_totals(void);
47 char *get_structure_title(void);
48 void vibrate(int natom,double *atomwt,double *x,double *y,double *z,double *charge);
49 void charge_dipole(int natom,double *x,double *y,double *z,double *atomwt,double *charge);
50 void set_connectivity(void);
51 // =====================================
52 void read_gmmxinp(char *paramfile,
53 FILE *infile,
54 int flags)
55 {
56 int i,ncount,nret;
57 int icount;
58 int ngood, nbparam, nbatom,nbcontact;
59 int batom = FALSE;
60 float logp;
61 //char *std_file;
62 char bad_atom[80],bad_contact[80],bad_param[80];
63 FILE *logfile;
64
65 // check filetype
66 ncount = 0;
67
68 /*
69 std_file = strdup(out_file);
70 strcpy(Savebox.fname, out_file);
71 */
72
73 // build filenames
74 strcpy(bad_atom,"_batom.sdf");
75 strcpy(bad_contact,"_bcontact.sdf");
76 strcpy(bad_param,"_bparam.sdf");
77
78 // open logfile
79 logfile = pcmlogfile;
80
81 files.nfiles = 0;
82 //
83 minim_control.field = MMFF94;
84 //strcpy(line,paramfile);
85 zero_data();
86 read_datafiles(paramfile);
87 minim_control.method = 3; // just do first deriv minimization
88 icount = 0;
89 //
90 ngood = nbparam = nbatom = nbcontact = 0;
91 i = 0;
92 while (1)
93 // for (i=1; i <= files.nfiles; i++)
94 {
95 ++i;
96 initialize();
97 nret = rd_sdf(infile);
98 if (-1 == nret) return ;
99 if (FALSE == nret)
100 {
101 fprintf(logfile,"Strnum: %d Cannot handle aromatic bonds. Try Kekule version.\n",i);
102 goto L_10;
103 natom = 0;
104 }
105
106 if (natom == 1)
107 {
108 fprintf(logfile,"Strnum: %d has only one atom\n",i);
109 goto L_10;
110 }
111 if (nret == TRUE)
112 {
113 // TJO
114 if (flags & DO_ADDH)
115 {
116 hdel(0);
117 hadd();
118 }
119 set_connectivity(); // count number of connections and total bond order for all atoms
120 type();
121
122 #ifdef KEG_DEBUG
123 minim_values.iprint = TRUE;
124 for(i=1; i<= natom; i++)
125 fprintf(pcmlogfile,"Atom: %d %d %ld\n",i,atom.type[i],atom.flags[i]);
126 #endif
127 nret = setup_calculation();
128
129 #ifdef KEG_DEBUG
130 print_energy_totals();
131 for(i=1; i<= natom; i++)
132 fprintf(pcmlogfile,"Atom: %d %d %ld\n",i,atom.type[i],atom.flags[i]);
133 if (minim_values.iprint == TRUE) exit(0);
134 minim_values.iprint = FALSE;
135 #endif
136 if (nret == TRUE)
137 {
138 minimize(natom,atom.use,atom.x,atom.y,atom.z);
139 // compute xlogp
140 if (flags & DO_XLOGP) {
141 if (VERBOSE) fprintf(pcmlogfile, " Evaluating XlogP\n");
142 logp = xlogp(natom,atom.atomnum,atom.iat,atom.bo,atom.flags,atom.type);
143 logp_calc.logp = logp;
144 }
145 // compute dipole moment
146 if (flags & DO_DIPOLE) {
147 if (VERBOSE) fprintf(pcmlogfile, " Evaluating dipole moment\n");
148 charge_dipole(natom,atom.x,atom.y,atom.z,atom.atomwt,atom.charge);
149 }
150 // compute vib
151 if (flags & DO_VIBRATION) {
152 if (VERBOSE) fprintf(pcmlogfile, " Evaluating vibrational parameters\n");
153 vibrate(natom,atom.atomwt,atom.x,atom.y,atom.z,atom.charge);
154 }
155 write_sdf(flags);
156 ++ngood;
157 ++files.nfiles;
158 #ifdef KEG_DEBUG
159 print_energy_totals();
160 #endif
161 } else
162 {
163 strcpy(Savebox.fname,bad_param);
164 files.append = TRUE;
165 //write_sdf(flags);
166 nbparam++;
167 //strcpy(Savebox.fname,std_file);
168 fprintf(logfile,"Strnum: %d CID: %s fails on missing parameters\n",files.nfiles,get_structure_title());
169 if (VERBOSE)
170 fprintf(stdout, "Missing parameters - no calc - %s\n",get_structure_title());
171 }
172 end_calculation();
173 if (VERBOSE) {
174 fprintf(stdout,
175 "Strnum: %d %s natom %d \n",files.nfiles,get_structure_title(), natom);
176 }
177 } else // read_sdf failed
178 {
179 fprintf(logfile,"Strnum: %d CID: %s fails on unsupported atom\n",files.nfiles,get_structure_title());
180 strcpy(Savebox.fname,bad_atom);
181 files.append = TRUE;
182 if (batom == FALSE)
183 {
184 batom = TRUE;
185 files.append = FALSE;
186 }
187 //write_sdf(flags);
188 nbatom++;
189 //strcpy(Savebox.fname,std_file);
190 }
191 L_10:
192 continue;
193 }
194 // fclose(infile);
195 fprintf(logfile,"Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",ngood,nbparam,nbatom,nbcontact);
196 if (VERBOSE) {
197 fprintf(stdout,
198 "Ngood: %d Nbparam: %d Nbatom: %d Nbcontact: %d \n",
199 ngood,nbparam,nbatom,nbcontact);
200 }
201
202 // fclose(logfile);
203 }
204 // ================================================
205 int close_contact(int mode)
206 {
207 int i, j;
208 double dx,dy,dz;
209
210 for (i=1; i < natom; i++)
211 {
212 for (j=i+1; j <= natom; j++)
213 {
214 dx = fabs(atom.x[i] - atom.x[j]);
215 dy = fabs(atom.y[i] - atom.y[j]);
216 dz = fabs(atom.z[i] - atom.z[j]);
217 if ( (dx + dy + dz) < 0.1)
218 {
219 return FALSE;
220 }
221 }
222 }
223 return TRUE;
224 }