ViewVC Help

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

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

#	Line 3 | Line 3
3
4		#include "pcwin.h"
5		#include "pcmod.h"
6	–	#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 27 | 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 35 | Line 40
40		void message_alert(char *, char *);
41		void hdel(int);
42		FILE * fopen_path ( char * , char * , char * ) ;
38	–	void read_schakal(int,int);
39	–	void write_schakal(void);
40	–	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 49 | 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 81 | Line 84
84		else
85		wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86		*/
84	–	wfile = stdout;
87
88		j = strlen(Struct_Title);
89		for (i=0; i < j; i++)
#	Line 152 | 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 230 | 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 274 | 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)
438	>	// =====================
439	>	void get_tag(char *line,char *tag)
440		{
441	<	int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
442	<	int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
443	<	int jji, jj1, jj2, jj3, jj4;
444	<	int has_Aromatic, xPlus,yPlus,zPlus, planar;
445	<	int icount,itemp[4];
446	<	float xtmp, ytmp, ztmp, dz;
447	<	char  c1[4],c2[4];
448	<	char  atomchar[3];
449	<	char  inputline[150];
450	<	FILE   *infile;
451	<
452	<	infile = fopen_path(Openbox.path,Openbox.fname,"r");
453	<
454	<	if (infile == NULL)
455	<	{
456	<	message_alert("Error opening Concord file","Concord Setup");
457	<	return FALSE;
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		}
307	–
308	–	if (isubred == 0)
309	–	ibotptr = 0;
310	–	else
311	–	{
312	–	ibotptr = natom;
313	–	substr.istract[substr.icurstr] = TRUE;
314	–	}
315	–
316	–	if (istruct > 1)
317	–	{
318	–	ncount = 0;
319	–	while (FetchRecord(infile,inputline))
320	–	{
321	–	if (strncasecmp(inputline,"$$$$",4) == 0)
322	–	{
323	–	ncount++;
324	–	if (ncount+1 == istruct)
325	–	goto L_1;
326	–	}
327	–	}
328	–	}
329	–	// start file read here
330	–	L_1:
331	–	got_title = FALSE;
332	–	xPlus = yPlus = zPlus = FALSE;
333	–	FetchRecord(infile,inputline);
334	–	sscanf(inputline,"SDF %s",Struct_Title);
335	–	got_title = TRUE;
336	–	/*     if (strlen(inputline) > 4)
337	–	{
338	–	iz = strlen(inputline);
339	–	if (iz > 60) iz = 59;
340	–	for (i=4; i < iz; i++)
341	–	Struct_Title[i-4] = inputline[i];
342	–	Struct_Title[i] = '\0';
343	–	got_title = TRUE;
344	–	} */
345	–	FetchRecord(infile,inputline); // blank line
346	–	FetchRecord(infile,inputline); // blank line
347	–
348	–	FetchRecord(infile,inputline); // natom and nbond
349	–	for (i=0; i <4; i++)
350	–	{
351	–	c1[i] = inputline[i];
352	–	c2[i] = inputline[i+3];
353	–	}
354	–	c1[3] = '\0';c2[3] = '\0';
355	–	niatom = atoi(c1);
356	–	nibond = atoi(c2);
357	–	if (niatom == 0)
358	–	return FALSE;
359	–
360	–	for (i=0; i < niatom; i++)
361	–	{
362	–	FetchRecord(infile,inputline);
363	–	sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
364	–
365	–	itype = 0;
366	–	if (xtmp != 0.0) xPlus= TRUE;
367	–	if (ytmp != 0.0) yPlus= TRUE;
368	–	if (ztmp != 0.0) zPlus= TRUE;
369	–
370	–	iz = strlen(atomchar);
371	–	if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1)  // nitro fix
372	–	{
373	–	if (atomchar[0] == 'N') itype = 8;
374	–	if (atomchar[0] == 'O') itype = 6;
375	–
376	–	if (junk1 != 0)
377	–	{
378	–	if (junk1 == 3 && itype == 8)
379	–	itype = 41;  // N+
380	–	if (junk1 == 5 && itype == 6)
381	–	itype = 42;  // O-
382	–	}
383	–	}
384	–
385	–	newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
386	–	if (newatom == -1)
387	–	{
388	–	printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
389	–	return FALSE;
390	–	}
391	–	if (isubred == 1)
392	–	{
393	–	atom[newatom].substr[0] = 0;
394	–	atom[newatom].substr[0] |= (1L << substr.icurstr);
395	–	}
396	–	}
397	–	has_Aromatic = FALSE;
398	–
399	–	planar = TRUE;
400	–	if (xPlus && yPlus && zPlus) planar = FALSE;
401	–
402	–	for (i=0; i < nibond; i++)
403	–	{
404	–	FetchRecord(infile,inputline);
405	–	for (j=0; j <4; j++)
406	–	{
407	–	c1[j] = inputline[j];
408	–	c2[j] = inputline[j+3];
409	–	}
410	–	c1[3] = '\0';c2[3] = '\0';
411	–	ia1 = atoi(c1);
412	–	ia2 = atoi(c2);
413	–	istereo = 0;
414	–	sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
415	–	if (ibond >= 4)
416	–	{
417	–	ibond = 1;
418	–	has_Aromatic = TRUE;
419	–	}
420	–	make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
421	–	if (istereo == 1 && planar)
422	–	{
423	–	atom[ia2+ibotptr].z += 0.7;
424	–	if (atom[ia2+ibotptr].atomnum == 1)
425	–	atom[ia2+ibotptr].type = 60;
426	–	}
427	–	if (istereo == 6 && planar)
428	–	{
429	–	atom[ia2+ibotptr].z -= 0.7;
430	–	if (atom[ia2+ibotptr].atomnum == 1)
431	–	atom[ia2+ibotptr].type = 60;
432	–	}
433	–	}
434	–	FetchRecord(infile,inputline);  // M END line
435	–	FetchRecord(infile,inputline);
436	–	if (got_title == FALSE)
437	–	strcpy(Struct_Title,inputline);
438	–	fclose(infile);
439	–	ncount = strlen(Struct_Title);
440	–	for (i=0; i < ncount; i++)
441	–	{
442	–	if (Struct_Title[i] == '\n')
443	–	{
444	–	Struct_Title[i] = '\0';
445	–	break;
446	–	}
447	–	}
448	–	//   search for quaternary centers - check if flat
449	–	dz = 0.0;
450	–	for (i=1; i <= natom; i++)
451	–	{
452	–	dz += atom[i].z;
453	–	jji = jj1 = jj2 = jj3 = jj4 = 0;
454	–	if (atom[i].type != 0)
455	–	{
456	–	for (j=0; j < 6; j++)
457	–	{
458	–	if (atom[i].iat[j] != 0)
459	–	jji++;
460	–	}
461	–	}
462	–	if (jji == 4)
463	–	{
464	–	if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
465	–	atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0)  // flat center
466	–	{
467	–	for (j=0; j < 4; j++)
468	–	{
469	–	itemp[j] = 0;
470	–	if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
471	–	if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
472	–	if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
473	–	if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
474	–	}
475	–	icount = 0;
476	–	if (jj1 == 1)
477	–	{
478	–	itemp[icount] = atom[i].iat[0];
479	–	icount++;
480	–	}
481	–	if (jj2 == 1)
482	–	{
483	–	itemp[icount] = atom[i].iat[1];
484	–	icount++;
485	–	}
486	–	if (jj3 == 1)
487	–	{
488	–	itemp[icount] = atom[i].iat[2];
489	–	icount++;
490	–	}
491	–	if (jj3 == 1)
492	–	{
493	–	itemp[icount] = atom[i].iat[3];
494	–	icount++;
495	–	}
496	–	if (icount >= 2)
497	–	{
498	–	atom[itemp[0]].z += 0.5;
499	–	atom[itemp[1]].z -= 0.5;
500	–	}
501	–	}
502	–	}
503	–	}
504	–	for (i=1; i < natom; i++)
505	–	{
506	–	for (j=i+1; j <= natom; j++)
507	–	{
508	–	xtmp = atom[i].x - atom[j].x;
509	–	ytmp = atom[i].y - atom[j].y;
510	–	ztmp = atom[i].z - atom[j].z;
511	–	if ( (xtmp + ytmp + ztmp) < 0.1)
512	–	{
513	–	atom[i].x += 0.1;
514	–	atom[i].y += 0.1;
515	–	atom[i].z += 0.1;
516	–	atom[j].x -= 0.1;
517	–	atom[j].y -= 0.1;
518	–	atom[j].z -= 0.1;
519	–	}
520	–	}
521	–	}
522	–
523	–	if (has_Aromatic == TRUE)
524	–	mopaco(1,natom);
525	–	return TRUE;
526	–	}
527	–	// ==================================================
528	–	void read_schakal(int isel,int isubred)
529	–	{
530	–	int i, itype,  newatom, iz;
531	–	float xtmp, ytmp, ztmp;
532	–	char  dummy[30];
533	–	char  atomchar[6];
534	–	char  inputline[150];
535	–	FILE   *infile;
536	–
537	–	infile = fopen_path(Openbox.path,Openbox.fname,"r");
538	–
539	–	if (infile == NULL)
540	–	{
541	–	message_alert("Error opening Schakal file","Schakal Setup");
542	–	return;
543	–	}
544	–	while ( FetchRecord(infile,inputline))
545	–	{
546	–	sscanf(inputline,"%s",dummy);
547	–	if (strcmp(dummy,"TITLE") == 0)
548	–	{
549	–	sscanf(inputline,"%s %s",dummy,Struct_Title);
550	–	} else if (strcmp(dummy,"ATOM") == 0)
551	–	{
552	–	sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
553	–	// strip off numbers
554	–	iz = strlen(atomchar);
555	–	for (i=0; i < iz; i++)
556	–	{
557	–	if (isdigit(atomchar[i]))
558	–	{
559	–	atomchar[i] = '\0';
560	–	break;
561	–	}
562	–	}
563	–	itype = 0;
564	–	newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
565	–	} else if (strcmp(dummy,"END") == 0)
566	–	{
567	–	break;
568	–	}
569	–	}
570	–	fclose(infile);
571	–	mopaco(1,natom);
572	–	//
462		}
574	–	// ==================================================
575	–	void write_schakal()
576	–	{
577	–	int i,j,lptest;
578	–	char atomchar[6];
579	–	FILE *wfile;
580	–
581	–	lptest = 0;
582	–	for( i = 1; i <= natom; i++ )
583	–	{
584	–	if( atom[i].mmx_type == 20 )
585	–	{
586	–	lptest = 1;
587	–	hdel( lptest );
588	–	break;
589	–	}
590	–	}
591	–
592	–	/*     now write the schakal file */
593	–	if (files.append)
594	–	wfile = fopen_path(Savebox.path,Savebox.fname,"a");
595	–	else
596	–	wfile = fopen_path(Savebox.path,Savebox.fname,"w");
597	–
598	–	j = strlen(Struct_Title);
599	–	for (i=0; i < j; i++)
600	–	{
601	–	if (Struct_Title[i] == '\n')
602	–	{
603	–	Struct_Title[i] = '\0';
604	–	break;
605	–	}
606	–	}
607	–	fprintf(wfile,"TITLE  %s\n",Struct_Title);
608	–	fprintf(wfile,"CELL \n");
609	–	for (i=1; i <= natom; i++)
610	–	{
611	–	sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
612	–
613	–	fprintf(wfile,"ATOM  %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
614	–	atom[i].y, atom[i].z);
615	–	}
616	–	fprintf(wfile,"END \n");
617	–	fclose(wfile);
618	–	}
619	–
620	–

Diff Legend

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