ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/freemol/trunk/src/mengine/src/read_sdf.c

#	Line 1 | Line 1
1	+	/* NOTICE: this source code file has been modified for use with FreeMOL */
2		#define EXTERN extern
3
4		#include "pcwin.h"
5		#include "pcmod.h"
5	–	#include "substr.h"
6		#include "atom_k.h"
7		#include "energies.h"
8	+	#include "fix.h"
9	+
10	+	#define TABULATOR_INCLUDE_IMPLEMENTATION
11	+	#include "tabulator.h"
12
13		EXTERN struct t_files {
14		int nfiles, append, batch, icurrent, ibatno;
#	Line 26 | Line 30
30		float etot,htot,stot,gtot,cptot;
31		} vibdata;
32
33	+	// ============================
34	+
35		int make_atom(int, float, float, float,char *);
36		void make_bond(int, int, int);
37		int read_sdf(int,int);
#	Line 34 | Line 40
40		void message_alert(char *, char *);
41		void hdel(int);
42		FILE * fopen_path ( char * , char * , char * ) ;
37	–	void read_schakal(int,int);
38	–	void write_schakal(void);
39	–	void mopaco(int,int);
43		int FetchRecord(FILE *, char *);
44		void avgleg(void);
45	<
45	>	void get_tag(char *,char *);
46	>	// ==================================
47		/*
48		* flags - indicates whether dipole, xlogp or vibrational calculations were done
49		* and if so write them to the output file. Look at the definitions in pcmod.h
#	Line 48 | Line 52
52		void write_sdf(int flags)
53		{
54		int i,j,lptest,nbond,junk;
55	<	FILE *wfile;
55	>	FILE *wfile = pcmoutfile;
56
57		lptest = 0;
58		for( i = 1; i <= natom; i++ )
#	Line 80 | Line 84
84		else
85		wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86		*/
83	–	wfile = stdout;
87
88		j = strlen(Struct_Title);
89		for (i=0; i < j; i++)
#	Line 151 | Line 154
154		fprintf(wfile,"$$$$\n");
155		//    fclose(wfile);
156		}
157	+
158		/* ===================================  */
159		// fast read - assume file is open and positioned
160		// read structure and down to end $$$$
161		// and return
162		int rd_sdf(FILE *rfile)
163		{
164	<	int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
165	<	int ncount,junk,junk1,got_title,istereo,iz;
164	>	int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom;
165	>	int ncount,junk,junk1,got_title,istereo,iz,nvalue;
166		int jji, jj1, jj2, jj3, jj4;
167	+	int n_row;
168		int has_Aromatic, xPlus,yPlus,zPlus, planar;
169		int icount,itemp[4];
170		int Ret_Val;
171		float xtmp, ytmp, ztmp, dz;
172	+	float fconst,min,max;
173		char  c1[4],c2[4];
174		char  atomchar[3];
175		char  inputline[150];
176	+	char  tag[30];
177	+	char ***tab;
178
179		Ret_Val = TRUE;
180		got_title = FALSE;
#	Line 229 | Line 237
237		atom[newatom].mmx_type = itype;
238		if (newatom == -1)
239		{
240	<	printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
240	>	fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
241		Ret_Val = FALSE;
242		}
243		}
#	Line 255 | Line 263
263		{
264		ibond = 1;
265		has_Aromatic = TRUE;
266	+	// TJO
267	+	Ret_Val = FALSE;
268		}
269		make_bond(ia1, ia2, ibond);
270		if (istereo == 1 && planar)
#	Line 271 | Line 281
281		}
282		}
283		// read to end of structure
284	+	// parse any tags here
285	+	// agreed tags
286	+	//  <FIXED_ATOMS>
287	+	//  <RESTRAINED_ATOMS>
288	+	//  <RESTRAINED_DISTANCES>
289	+	//  <RESTRAINED_ANGLES>
290	+	//  <RESTRAINED_DIHEDRALS>
291	+
292		while (FetchRecord(rfile,inputline))
293		{
294		if (strncasecmp(inputline,"$$$$",4) == 0)
295		{
296		return Ret_Val;
297	+	} else if (inputline[0] == '>')
298	+	{
299	+	get_tag(inputline,tag);
300	+	if (strcasecmp(tag,"fixed_atoms") == 0)
301	+	{
302	+	tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
303	+	if (tab)
304	+	{
305	+	fprintf(pcmlogfile,"got table\n");
306	+	n_row = tabulator_height(tab);
307	+	for (i=0; i < n_row; i++)
308	+	{
309	+	ia1 = atoi(tab[i][0]);
310	+	if (ia1 > 0 && ia1 <= natom)
311	+	{
312	+	fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
313	+	fx_atom.natom_fix++;
314	+	}
315	+	}
316	+	}
317	+	} else if (strcasecmp(tag,"restrained_atoms") == 0)
318	+	{
319	+	tab = tabulator_new_from_file_using_header(rfile,"ATOM1 MAX F_CONST X Y Z",0);
320	+	if (tab)
321	+	{
322	+	n_row = tabulator_height(tab);
323	+	for (i=0; i < n_row; i++)
324	+	{
325	+	ia1 = atoi(tab[i][0]);
326	+	max = atof(tab[i][1]);
327	+	fconst = atof(tab[i][2]);
328	+	if (ia1 > 0 && ia1 <= natom)
329	+	{
330	+	restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
331	+	restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
332	+	restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
333	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
334	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
335	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
336	+	if (strlen(tab[i][3]) > 0)  // x postion
337	+	{
338	+	xtmp = atof(tab[i][3]);
339	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
340	+	}
341	+	if (strlen(tab[i][4]) > 0)  // y postion
342	+	{
343	+	xtmp = atof(tab[i][4]);
344	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
345	+	}
346	+	if (strlen(tab[i][5]) > 0)  // y postion
347	+	{
348	+	xtmp = atof(tab[i][3]);
349	+	restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
350	+	}
351	+	restrain_atom.natom_restrain++;
352	+	}
353	+	}
354	+	}
355	+	} else if (strcasecmp(tag,"restrained_distances") == 0)
356	+	{
357	+	tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
358	+	if (tab)
359	+	{
360	+	n_row = tabulator_height(tab);
361	+	for (i=0; i < n_row; i++)
362	+	{
363	+	ia1 = atoi(tab[i][0]);
364	+	ia2 = atoi(tab[i][1]);
365	+	min = atof(tab[i][2]);
366	+	max = atof(tab[i][3]);
367	+	fconst = atof(tab[i][4]);
368	+	if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
369	+	{
370	+	fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
371	+	fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
372	+	fx_dist.fdconst[fx_dist.ndfix] = fconst;
373	+	fx_dist.min_dist[fx_dist.ndfix] = min;
374	+	fx_dist.max_dist[fx_dist.ndfix] = max;
375	+	fx_dist.ndfix++;
376	+	}
377	+	}
378	+	}
379	+	} else if (strcasecmp(tag,"restrained_angles") == 0)
380	+	{
381	+	tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
382	+	if (tab)
383	+	{
384	+	n_row = tabulator_height(tab);
385	+	for (i=0; i < n_row; i++)
386	+	{
387	+	ia1 = atoi(tab[i][0]);
388	+	ia2 = atoi(tab[i][1]);
389	+	ia3 = atoi(tab[i][2]);
390	+	min = atof(tab[i][3]);
391	+	max = atof(tab[i][4]);
392	+	fconst = atof(tab[i][5]);
393	+	if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
394	+	{
395	+	fx_angle.kafix[fx_angle.nafix][0] = ia1;
396	+	fx_angle.kafix[fx_angle.nafix][1] = ia2;
397	+	fx_angle.kafix[fx_angle.nafix][2] = ia3;
398	+	fx_angle.faconst[fx_angle.nafix] = fconst;
399	+	fx_angle.min_ang[fx_angle.nafix] = min;
400	+	fx_angle.max_ang[fx_angle.nafix] = max;
401	+	fx_angle.nafix++;
402	+	}
403	+	}
404	+	}
405	+	} else if (strcasecmp(tag,"restrained_dihedrals") == 0)
406	+	{
407	+	tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
408	+	if (tab)
409	+	{
410	+	n_row = tabulator_height(tab);
411	+	for (i=0; i < n_row; i++)
412	+	{
413	+	ia1 = atoi(tab[i][0]);
414	+	ia2 = atoi(tab[i][1]);
415	+	ia3 = atoi(tab[i][2]);
416	+	ia4 = atoi(tab[i][3]);
417	+	min = atof(tab[i][4]);
418	+	max = atof(tab[i][5]);
419	+	fconst = atof(tab[i][6]);
420	+	if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
421	+	{
422	+	fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
423	+	fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
424	+	fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
425	+	fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
426	+	fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
427	+	fx_torsion.min_tor[fx_torsion.ntfix] = min;
428	+	fx_torsion.max_tor[fx_torsion.ntfix] = max;
429	+	fx_torsion.ntfix++;
430	+	}
431	+	}
432	+	}
433	+	}
434		}
435		}
436		return Ret_Val;
437		}
438	<	/* ===================================  */
439	<	int read_sdf(int istruct, int isubred)
285	<	{
286	<	int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
287	<	int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
288	<	int jji, jj1, jj2, jj3, jj4;
289	<	int has_Aromatic, xPlus,yPlus,zPlus, planar;
290	<	int icount,itemp[4];
291	<	float xtmp, ytmp, ztmp, dz;
292	<	char  c1[4],c2[4];
293	<	char  atomchar[3];
294	<	char  inputline[150];
295	<	FILE   *infile;
296	<
297	<	infile = fopen_path(Openbox.path,Openbox.fname,"r");
298	<
299	<	if (infile == NULL)
300	<	{
301	<	message_alert("Error opening Concord file","Concord Setup");
302	<	return FALSE;
303	<	}
304	<
305	<	if (isubred == 0)
306	<	ibotptr = 0;
307	<	else
308	<	{
309	<	ibotptr = natom;
310	<	substr.istract[substr.icurstr] = TRUE;
311	<	}
312	<
313	<	if (istruct > 1)
314	<	{
315	<	ncount = 0;
316	<	while (FetchRecord(infile,inputline))
317	<	{
318	<	if (strncasecmp(inputline,"$$$$",4) == 0)
319	<	{
320	<	ncount++;
321	<	if (ncount+1 == istruct)
322	<	goto L_1;
323	<	}
324	<	}
325	<	}
326	<	// start file read here
327	<	L_1:
328	<	got_title = FALSE;
329	<	xPlus = yPlus = zPlus = FALSE;
330	<	FetchRecord(infile,inputline);
331	<	sscanf(inputline,"SDF %s",Struct_Title);
332	<	got_title = TRUE;
333	<	/*     if (strlen(inputline) > 4)
334	<	{
335	<	iz = strlen(inputline);
336	<	if (iz > 60) iz = 59;
337	<	for (i=4; i < iz; i++)
338	<	Struct_Title[i-4] = inputline[i];
339	<	Struct_Title[i] = '\0';
340	<	got_title = TRUE;
341	<	} */
342	<	FetchRecord(infile,inputline); // blank line
343	<	FetchRecord(infile,inputline); // blank line
344	<
345	<	FetchRecord(infile,inputline); // natom and nbond
346	<	for (i=0; i <4; i++)
347	<	{
348	<	c1[i] = inputline[i];
349	<	c2[i] = inputline[i+3];
350	<	}
351	<	c1[3] = '\0';c2[3] = '\0';
352	<	niatom = atoi(c1);
353	<	nibond = atoi(c2);
354	<	if (niatom == 0)
355	<	return FALSE;
356	<
357	<	for (i=0; i < niatom; i++)
358	<	{
359	<	FetchRecord(infile,inputline);
360	<	sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
361	<
362	<	itype = 0;
363	<	if (xtmp != 0.0) xPlus= TRUE;
364	<	if (ytmp != 0.0) yPlus= TRUE;
365	<	if (ztmp != 0.0) zPlus= TRUE;
366	<
367	<	iz = strlen(atomchar);
368	<	if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1)  // nitro fix
369	<	{
370	<	if (atomchar[0] == 'N') itype = 8;
371	<	if (atomchar[0] == 'O') itype = 6;
372	<
373	<	if (junk1 != 0)
374	<	{
375	<	if (junk1 == 3 && itype == 8)
376	<	itype = 41;  // N+
377	<	if (junk1 == 5 && itype == 6)
378	<	itype = 42;  // O-
379	<	}
380	<	}
381	<
382	<	newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
383	<	if (newatom == -1)
384	<	{
385	<	printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
386	<	return FALSE;
387	<	}
388	<	if (isubred == 1)
389	<	{
390	<	atom[newatom].substr[0] = 0;
391	<	atom[newatom].substr[0] |= (1L << substr.icurstr);
392	<	}
393	<	}
394	<	has_Aromatic = FALSE;
395	<
396	<	planar = TRUE;
397	<	if (xPlus && yPlus && zPlus) planar = FALSE;
398	<
399	<	for (i=0; i < nibond; i++)
400	<	{
401	<	FetchRecord(infile,inputline);
402	<	for (j=0; j <4; j++)
403	<	{
404	<	c1[j] = inputline[j];
405	<	c2[j] = inputline[j+3];
406	<	}
407	<	c1[3] = '\0';c2[3] = '\0';
408	<	ia1 = atoi(c1);
409	<	ia2 = atoi(c2);
410	<	istereo = 0;
411	<	sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
412	<	if (ibond >= 4)
413	<	{
414	<	ibond = 1;
415	<	has_Aromatic = TRUE;
416	<	}
417	<	make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
418	<	if (istereo == 1 && planar)
419	<	{
420	<	atom[ia2+ibotptr].z += 0.7;
421	<	if (atom[ia2+ibotptr].atomnum == 1)
422	<	atom[ia2+ibotptr].type = 60;
423	<	}
424	<	if (istereo == 6 && planar)
425	<	{
426	<	atom[ia2+ibotptr].z -= 0.7;
427	<	if (atom[ia2+ibotptr].atomnum == 1)
428	<	atom[ia2+ibotptr].type = 60;
429	<	}
430	<	}
431	<	FetchRecord(infile,inputline);  // M END line
432	<	FetchRecord(infile,inputline);
433	<	if (got_title == FALSE)
434	<	strcpy(Struct_Title,inputline);
435	<	fclose(infile);
436	<	ncount = strlen(Struct_Title);
437	<	for (i=0; i < ncount; i++)
438	<	{
439	<	if (Struct_Title[i] == '\n')
440	<	{
441	<	Struct_Title[i] = '\0';
442	<	break;
443	<	}
444	<	}
445	<	//   search for quaternary centers - check if flat
446	<	dz = 0.0;
447	<	for (i=1; i <= natom; i++)
448	<	{
449	<	dz += atom[i].z;
450	<	jji = jj1 = jj2 = jj3 = jj4 = 0;
451	<	if (atom[i].type != 0)
452	<	{
453	<	for (j=0; j < 6; j++)
454	<	{
455	<	if (atom[i].iat[j] != 0)
456	<	jji++;
457	<	}
458	<	}
459	<	if (jji == 4)
460	<	{
461	<	if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
462	<	atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0)  // flat center
463	<	{
464	<	for (j=0; j < 4; j++)
465	<	{
466	<	itemp[j] = 0;
467	<	if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
468	<	if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
469	<	if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
470	<	if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
471	<	}
472	<	icount = 0;
473	<	if (jj1 == 1)
474	<	{
475	<	itemp[icount] = atom[i].iat[0];
476	<	icount++;
477	<	}
478	<	if (jj2 == 1)
479	<	{
480	<	itemp[icount] = atom[i].iat[1];
481	<	icount++;
482	<	}
483	<	if (jj3 == 1)
484	<	{
485	<	itemp[icount] = atom[i].iat[2];
486	<	icount++;
487	<	}
488	<	if (jj3 == 1)
489	<	{
490	<	itemp[icount] = atom[i].iat[3];
491	<	icount++;
492	<	}
493	<	if (icount >= 2)
494	<	{
495	<	atom[itemp[0]].z += 0.5;
496	<	atom[itemp[1]].z -= 0.5;
497	<	}
498	<	}
499	<	}
500	<	}
501	<	for (i=1; i < natom; i++)
502	<	{
503	<	for (j=i+1; j <= natom; j++)
504	<	{
505	<	xtmp = atom[i].x - atom[j].x;
506	<	ytmp = atom[i].y - atom[j].y;
507	<	ztmp = atom[i].z - atom[j].z;
508	<	if ( (xtmp + ytmp + ztmp) < 0.1)
509	<	{
510	<	atom[i].x += 0.1;
511	<	atom[i].y += 0.1;
512	<	atom[i].z += 0.1;
513	<	atom[j].x -= 0.1;
514	<	atom[j].y -= 0.1;
515	<	atom[j].z -= 0.1;
516	<	}
517	<	}
518	<	}
519	<
520	<	if (has_Aromatic == TRUE)
521	<	mopaco(1,natom);
522	<	return TRUE;
523	<	}
524	<	// ==================================================
525	<	void read_schakal(int isel,int isubred)
438	>	// =====================
439	>	void get_tag(char *line,char *tag)
440		{
441	<	int i, itype,  newatom, iz;
442	<	float xtmp, ytmp, ztmp;
443	<	char  dummy[30];
444	<	char  atomchar[6];
445	<	char  inputline[150];
446	<	FILE   *infile;
447	<
448	<	infile = fopen_path(Openbox.path,Openbox.fname,"r");
449	<
450	<	if (infile == NULL)
451	<	{
452	<	message_alert("Error opening Schakal file","Schakal Setup");
453	<	return;
441	>	int i,j,icount;
442	>	icount = 0;
443	>	strcpy(tag,"");
444	>	for (i=0; i < strlen(line); i++)
445	>	{
446	>	if (line[i] == '<')
447	>	{
448	>	for (j=i+1; j < strlen(line); j++)
449	>	{
450	>	if (line[j] == '>')
451	>	{
452	>	tag[icount] = '\0';
453	>	return;
454	>	} else
455	>	{
456	>	tag[icount] = line[j];
457	>	icount++;
458	>	}
459	>	}
460	>	}
461		}
541	–	while ( FetchRecord(infile,inputline))
542	–	{
543	–	sscanf(inputline,"%s",dummy);
544	–	if (strcmp(dummy,"TITLE") == 0)
545	–	{
546	–	sscanf(inputline,"%s %s",dummy,Struct_Title);
547	–	} else if (strcmp(dummy,"ATOM") == 0)
548	–	{
549	–	sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
550	–	// strip off numbers
551	–	iz = strlen(atomchar);
552	–	for (i=0; i < iz; i++)
553	–	{
554	–	if (isdigit(atomchar[i]))
555	–	{
556	–	atomchar[i] = '\0';
557	–	break;
558	–	}
559	–	}
560	–	itype = 0;
561	–	newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
562	–	} else if (strcmp(dummy,"END") == 0)
563	–	{
564	–	break;
565	–	}
566	–	}
567	–	fclose(infile);
568	–	mopaco(1,natom);
569	–	//
462		}
571	–	// ==================================================
572	–	void write_schakal()
573	–	{
574	–	int i,j,lptest;
575	–	char atomchar[6];
576	–	FILE *wfile;
577	–
578	–	lptest = 0;
579	–	for( i = 1; i <= natom; i++ )
580	–	{
581	–	if( atom[i].mmx_type == 20 )
582	–	{
583	–	lptest = 1;
584	–	hdel( lptest );
585	–	break;
586	–	}
587	–	}
588	–
589	–	/*     now write the schakal file */
590	–	if (files.append)
591	–	wfile = fopen_path(Savebox.path,Savebox.fname,"a");
592	–	else
593	–	wfile = fopen_path(Savebox.path,Savebox.fname,"w");
594	–
595	–	j = strlen(Struct_Title);
596	–	for (i=0; i < j; i++)
597	–	{
598	–	if (Struct_Title[i] == '\n')
599	–	{
600	–	Struct_Title[i] = '\0';
601	–	break;
602	–	}
603	–	}
604	–	fprintf(wfile,"TITLE  %s\n",Struct_Title);
605	–	fprintf(wfile,"CELL \n");
606	–	for (i=1; i <= natom; i++)
607	–	{
608	–	sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
609	–
610	–	fprintf(wfile,"ATOM  %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
611	–	atom[i].y, atom[i].z);
612	–	}
613	–	fprintf(wfile,"END \n");
614	–	fclose(wfile);
615	–	}
616	–
617	–

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines