#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define NFORWARD 4
#define DIM 3
#define MAXBUF 5000
#define MAXCHAR_LINE  10000
#define READ(a,b)   fread(&(b),sizeof(b),1,(a))

inline int strindex(const char *a, const char *b, int isleft)
{
int la=(int)strlen(a);
int lb=(int)strlen(b);
  if (lb==0) return -1;
  if (lb>la) return -1;
int i;

 for (i=isleft?0:la-lb;isleft?i<=la-lb:i>=0;i+=isleft?1:-1)
  if (strncmp(a+i,b,lb)==0)
    return i;
 return -1;
}

void get_line (FILE *fp, char *buffer)
{
	
	
	int i=0;
	char gc;
	
	
	
	
	if (feof(fp)==0) READ(fp,gc);
	while ((gc!=10) && (feof(fp)==0)) 
	{
        if (gc!=13) { buffer[i]=gc;i++;}
        READ(fp,gc);
		if (i>MAXCHAR_LINE) {i=0;break;}       
	}
	
	buffer[i]=0;
	
	
	
}

#define MAXBUFFER 5000

void trcspace(char *tstr)
{
	char   tempstr[MAXBUFFER],temp[MAXBUFFER];
	int ti,tend,tspcnt;
	int i;
	
	
	if (strlen(tstr)==0) {
		tempstr[0]='|';
		tempstr[1]='n';
		tempstr[2]='o';
		tempstr[3]='n';
		tempstr[4]='e';
		i=5;
		
    }
	else
	{
		i=0;
		tend=(int)strlen(tstr);
		tspcnt=0;
		for( ti=0;ti< tend;ti++)  if (tstr[ti]==' ')  tspcnt++;
		if (tspcnt==tend) 
		{
			tempstr[0]='|';
			tempstr[1]='n';
			tempstr[2]='o';
			tempstr[3]='n';
			tempstr[4]='e';
			i=5;
			
		}
		else
		{
			int j;
			for (j=1;j<=tend;j++) 
				
			{temp[j]=tstr[j-1];
				
				if (temp[j]==9) temp[j]=' ';
			}
			temp[0]=' ';
			
			temp[tend+1]=0;     
			for( ti=0;ti<tend+1;ti++) 
				if (temp[ti]!=' ') 
				{
					if (temp[ti-1]==' ') { tempstr[i]='|';i++;}
					tempstr[i]=temp[ti];
					i++;
					
				}
			
			
			
			
		}
		
		
	}
	
	tempstr[i]='|';
	i++;
	tempstr[i]=0;
	
	for (ti=0;ti<=i;ti++) tstr[ti]=tempstr[ti]; 
	
	
}

#undef MAXBUFFER




int parcnt(char *pstr)
{
	int  ti,pj;
	
	pj=0;
	for(ti=0;ti<(int)strlen(pstr);ti++)  if (pstr[ti]=='|')  pj++;
	pj--;
	return pj;
	
	
	
	
}


int iseq(float a, float b, float tole)
{
  if (a<0.) a=-a;
  if (b<0.) b=-b;
  if (a==0.0&&b==0.0) return 1;
  if (a>b)
   {
     float c=b;
	b=a; a=c;
   }
   float d=(b-a)/b;
  if (d<tole) return 1; else return 0;
}



void parstr(char *pstr, int pnum, char *tempstr)
{
	int pai,pj;
	
	
	
	if (pnum>parcnt(pstr)) 
	{tempstr[0]=0;exit;}
	else
	{
		int i=0;
		
		
		pj=0;
		for( pai=0;pai<(int)strlen(pstr);pai++) 
			if (pstr[pai]=='|') { pj++;}
			else
			{
				if (pj==pnum) { tempstr[i]=pstr[pai];i++;}
				
			}
		
		tempstr[i]=0;
	}
	
	
}




int main (int agrc, char * argv[])
{
	
	int i,k,j,l,m,n;
	char buffer[MAXBUF],buf[MAXBUF];
#define MAXINTTYP 256
#define NAMINO 28
	int nparam=0;
	float vdwParam[MAXINTTYP][2];
	char aanm[NAMINO][11]={
		"ALA","GLY","PRO","DPR"
		,"LEU","LYS","ARG","MET"
		,"GLU","GLUH","ASP","ASPH"
		,"GLN","ASN","VAL","ILE"
		,"PHE","TYR","TRP","HIS"
		,"HISH","CYS","CYSH","SER"
		,"THR","ACE","NMA","HISE"
	};

#define NLIP 2
  char lipnm[NLIP][11]={
    "POPE","POPC"   

   };
	
	FILE *fp,*ftmp,*fprm;
	
	fp=fopen(argv[1],"r");
	int ok=0;
	
	
	while ((feof(fp)==0)&&(ok==0))
	{
		get_line(fp,buffer);
		trcspace (buffer);
		if (parcnt(buffer)==2)
		{
			parstr(buffer,2,buf);
			if (strcmp(buf,"!NATOM")==0)
			{ok=1;break;}
		}
		
	}
	
	if (ok==0)
	{
		printf( "Wrong PSF file!\n");
		exit (1);
	}
	
	parstr (buffer,1,buf);
	int psfnatom;
	sscanf (buf,"%d",&psfnatom);
	
	char **psfnewnm,**psfsegnm, **psfresnm, **psfatomnm, **psfatomtype;
	int *psfresid,*psfispair;
        int *islip;	
        int nlipatom=0;	
	/*0 not in the pair list 1 yes*/
	//loading all atom info from psf
	psfnewnm =malloc (psfnatom * sizeof (char *));
	psfsegnm = malloc (psfnatom * sizeof (char *));
	psfresnm = malloc (psfnatom * sizeof (char *));
	psfatomnm = malloc (psfnatom * sizeof (char *));
	psfatomtype = malloc (psfnatom * sizeof (char *));
	psfresid = malloc (psfnatom * sizeof (int));
        islip = malloc(psfnatom*sizeof(int));
       // for (i=0;i<psfnatom;islip[i++]=-1);
	int read_count;
	float x,y;
	
	for (i=0;i<psfnatom;i++)
	{
        psfnewnm[i] = malloc (11*sizeof (char));
	psfsegnm[i] = malloc (11*sizeof (char));
        psfresnm[i] = malloc (11*sizeof (char));
        psfatomnm[i] = malloc (11*sizeof (char));
        psfatomtype[i] = malloc (11*sizeof (char));
		psfnewnm[i][0]=0;
        islip[i]=-1;
	}
	
	int nproatom=0;
	
	for (i=0;i<psfnatom;i++)
	{
		get_line(fp,buffer);
		read_count=sscanf(buffer, "%d %s %d %s %s %s %f %f",
						  &j, psfsegnm[i], &(psfresid[i]),
						  psfresnm[i], psfatomnm[i], psfatomtype[i], &x, &y);
		
		if (read_count!=8)
		{
			printf( "Incorrect line for atom!\n");
			exit (1);
			
		}
		
		
		for (m=0;m<NAMINO;m++)
			if (strcmp(psfresnm[i],aanm[m])==0)
			{
				nproatom++;break;
			}
		
                for (m=0;m<NLIP;m++)
                        if (strcmp(psfresnm[i],lipnm[m])==0)
                        {
				islip[i]= nlipatom++;
                                break;
                        }

		
	}


	printf ("Info : %d protein atoms and %d lipid atoms identified\n", nproatom, nlipatom);
	
	int *idx2res;
	int npep;
	int nres=nproatom;
	
	int *HD21,*HD22,*ND2,*HD1,*ND1,*OD1,*OD2,*OG1, *CD1,*CD2, *OE, *NE, *AUX,*OG,*CG,*CG1,*CG2,*H, *CA,  *C, *NH, *CB, *O;
	char **resnm;
	
	idx2res= malloc ((nres+1)*sizeof(int));
	CA=malloc ((nres+1)*sizeof(int));
	C=malloc ((nres+1)*sizeof(int));
	NH=malloc ((nres+1)*sizeof(int));
	CB=malloc ((nres+1)*sizeof(int));
	O=malloc ((nres+1)*sizeof(int));
	H=malloc ((nres+1)*sizeof(int));
	CG1=malloc ((nres+1)*sizeof(int));
	CG2=malloc ((nres+1)*sizeof(int));
	CG=malloc ((nres+1)*sizeof(int));
	OG=malloc ((nres+1)*sizeof(int));
	CD1=malloc ((nres+1)*sizeof(int));
	CD2=malloc ((nres+1)*sizeof(int));
	OE=malloc ((nres+1)*sizeof(int));
	NE=malloc ((nres+1)*sizeof(int));
	AUX=malloc ((nres+1)*sizeof(int));
	OG1=malloc ((nres+1)*sizeof(int));
	OD1=malloc ((nres+1)*sizeof(int));
	OD2=malloc ((nres+1)*sizeof(int));
	HD1=malloc ((nres+1)*sizeof(int));
	ND1=malloc ((nres+1)*sizeof(int));
	ND2=malloc ((nres+1)*sizeof(int));
	HD21=malloc ((nres+1)*sizeof(int));
	HD22=malloc ((nres+1)*sizeof(int));
	resnm=malloc ((nres+1)*sizeof(char*));
	for (i=0;i<nres+1;i++)
	{
		CA[i]=C[i]=NH[i]=CB[i]=O[i]=
		H[i]=CG1[i]=CG2[i]=CG[i]=
		OG[i]=CD1[i]=CD2[i]=OE[i]=
		NE[i]=AUX[i]=OG1[i]=OD1[i]=
		OD2[i]=HD1[i]=ND1[i]=ND2[i]=
		HD21[i]=HD22[i]=idx2res[i]=-1;
		resnm[i]=malloc (11*sizeof(char));
		resnm[i][0]='#';
		resnm[i][1]=0;
	}
	
	
	buf[0]=0;
	i=0;
	nres=-1; /*for res*/
	npep=-1; /*for seg*/
	l=-1; /*for resid*/
	ok=0;
	while ((i<psfnatom))
	{
		if (strcmp(psfsegnm[i],buf)!=0)
		/*different segment*/
		{
			strcpy(buf,psfsegnm[i]);
			ok=0;
			for (m=0;m<NAMINO;m++)
				if (strcmp(psfresnm[i],aanm[m])==0)
				{
					ok=1;break;
				}

                        printf("Info : new segment %s start t %d\n", buf, i);
                        if (ok) printf("Info : This is a protein segment.\n");

			if (ok==0) /*non protein*/ {i++;continue;}
			
			/*new protein segment (also new res)*/
			
			l=psfresid[i];
			npep++;
			nres+=5;/*a blank residue between two pep segments*/
			idx2res[npep]=nres;
			strcpy(resnm[nres],psfresnm[i]);
			
			
			
			
		}	
		
		
		
		if (ok==0)  {i++;continue;}
		
		if (l!=psfresid[i]) {nres++;l=psfresid[i];strcpy(resnm[nres],psfresnm[i]);}
		
		if (strcmp(psfatomnm[i],"CA")==0)  CA[nres]=i;
        if (strcmp(psfatomnm[i],"C")==0)  C[nres]=i;
        if (strcmp(psfatomnm[i],"O")==0)  O[nres]=i;
        if (strcmp(psfatomnm[i],"N")==0)  NH[nres]=i;
        if (strcmp(psfatomnm[i],"H")==0)  H[nres]=i;
        if (strcmp(psfatomnm[i],"CB")==0)  CB[nres]=i;
        if (strcmp(psfatomnm[i],"CG")==0)  CG[nres]=i;
        if (strcmp(psfatomnm[i],"CG1")==0)  CG1[nres]=i;
        if (strcmp(psfatomnm[i],"CG2")==0)  CG2[nres]=i;
        if (strcmp(psfatomnm[i],"OG")==0)  OG[nres]=i;
        if (strcmp(psfatomnm[i],"OE1")==0)  OE[nres]=i;
        if (strcmp(psfatomnm[i],"NE2")==0)  NE[nres]=i;
        if (strcmp(psfatomnm[i],"CD1")==0)  CD1[nres]=i;
        if (strcmp(psfatomnm[i],"CD2")==0)  CD2[nres]=i;
        if (strcmp(psfatomnm[i],"AUXL")==0)  AUX[nres]=i;
        if (strcmp(psfatomnm[i],"OG1")==0)  OG1[nres]=i;
        if (strcmp(psfatomnm[i],"OD1")==0)  OD1[nres]=i;
        if (strcmp(psfatomnm[i],"OD2")==0)  OD2[nres]=i;
        if (strcmp(psfatomnm[i],"HD1")==0)  HD1[nres]=i;
        if (strcmp(psfatomnm[i],"ND1")==0)  ND1[nres]=i;
        if (strcmp(psfatomnm[i],"ND2")==0)  ND2[nres]=i;
        if (strcmp(psfatomnm[i],"HD21")==0)  HD21[nres]=i;
        if (strcmp(psfatomnm[i],"HD22")==0)  HD22[nres]=i;
		
		
		
		
		
		
		
		
		i++;	
		
	}
	
	
	if ((npep==-1)||(nres==-1))
	{
		printf("No protein atom found!\n");
		exit(1);
	}
	
	nres+=2;
	npep++;
	idx2res[npep]=nres;

        //for(m=0;m<npep;m++) printf("Info: protein segment %d  residue id range [%d - %d]\n"
        //                           ,m, idx2res[m], idx2res[m+1]);

	
	int  totres=nres;
	int npair=0;
	ftmp = fopen("temp.xvg","w");
	/*generate temp file for a pair list*/
        /*
	for (i=0;i<=npep;i++) fprintf(ftmp,"%d %d\n",i,idx2res[i]);
	for (i=0;i<nres;i++)
		fprintf(ftmp,"%s CA:%d C:%d O:%d N:%d H:%d CB:%d CG1:%d CG2:%d CG:%d OG:%d OG1:%d OE:%d NE:%d CD1:%d CD2:%d OD1:%d OD2:%d HD1:%d ND1:%d ND2:%d\n", 
				resnm[i],CA[i], C[i], O[i], NH[i], H[i],CB[i],CG1[i],CG2[i],CG[i],OG[i],OG1[i],OE[i],NE[i],CD1[i],CD2[i],OD1[i],OD2[i],HD1[i],ND1[i],ND2[i]);
	
	
	fprintf (ftmp,"[ pairs ]\n");
        */
	for (m=0;m<npep;m++)
	{
		/*segment m start from res idx2res[m]<= i < idx2res[m+1]*/
		nres = idx2res[m+1];
		
		for (i=idx2res[m];i<nres;i++)
		{
			
			
			
			
			//cai-->cbi+4   0.671  0.47
			
			if (i<=nres-4)
			{
				
				if ((CA[i]>=0) && (CB[i+4]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  2.89314E-02 3.11858E-04  ;  0.671  0.47 nm\n", CA[i]+1,  CB[i+4]+1);
					
				}
				
				
				if ((CA[i]>=0) && (CA[i+4]>=0)&&(strcmp(resnm[i+4],"GLY")==0))
				{ 
					npair++;fprintf(ftmp,"%d    %d  1  4.19375E-02  6.55273E-04  ;  0.671  0.50 GLY nm\n", CA[i]+1,  CA[i+4]+1);
					
				}
				
				
				
				
			}
			
			
			
			
			
			
			//reinforce Oi-Ni+3 HB for turn
			
			if (i<=nres-3)
			{
				if ((strcmp(resnm[i+3],"PRO")!=0)&&(strcmp(resnm[i+3],"DPR")!=0))
					if ((O[i]>=0) && (NH[i+3]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  1.12369E-02     2.14740E-06 ;  Oi-Ni+3 \n", O[i]+1,  NH[i+3]+1);
						
					}
				
				
				
			}
			
			
			
			
			if (i<=nres-4)
			{
				if ((strcmp(resnm[i+4],"PRO")!=0)&&(strcmp(resnm[i+4],"DPR")!=0))
					if ((O[i]>=0) && (NH[i+4]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  1.12369E-02     2.14740E-06 ;  Oi-Ni+4 \n", O[i]+1,  NH[i+4]+1);
						
					}
				
				
				
			}
			
			
			
			
			
			
			// correct for gamma particles
			
			if ((OG[i]>=0) && (H[i]>=0))
			{ 
				npair++;fprintf(ftmp,"%d    %d  1  1.13380E-03  1.28550E-07 ; 2.5 0.22 nm OG-H \n", OG[i]+1,  H[i]+1);
				
			}
			if ((OG[i]>=0) && (O[i]>=0))
			{ 
				npair++;fprintf(ftmp,"%d    %d  1 7.74841E-04   1.50095E-06 ; 2.5  0.27 nm OG-Oi  rep\n", OG[i]+1,  O[i]+1);
				
			}
			if ((OG[i]>=0) && (NH[i+1]>=0))
			{ 
				npair++;fprintf(ftmp,"%d    %d  1  3.83970E-03  4.12285E-06 ;   OG-Ni+1 \n", OG[i]+1,  NH[i+1]+1);
				
			}
			
			
			if ((OG[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.13380E-03  1.28550E-07  ;  Og-Hi+1\n", OG[i]+1,  H[i+1]+1);
				
			}
			
			
			
			if (i>idx2res[m])
			{
				
				if ((OG[i]>=0) && (O[i-1]>=0))
				{ 
					npair++;fprintf(ftmp,"%d    %d  1  3.87420E-03   1.50095E-06 ;  Og-Oi-1 \n", OG[i]+1,  O[i-1]+1);
					
				}
				
				
				
			}
			
			
			
			// og1
			if ((OG1[i]>=0) && (H[i]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  2.26760E-04  1.28550E-07  ; OG1-Hi rep\n", OG1[i]+1,  H[i]+1);
				
			}
			if ((OG1[i]>=0) && (O[i]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  7.74841E-04   1.50095E-06 ; OG1-Oi rep\n", OG1[i]+1,  O[i]+1);
				
			}
			if ((OG1[i]>=0) && (NH[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  3.83970E-03  4.12285E-06 ;  OG1-Ni+1 \n", OG1[i]+1,  NH[i+1]+1);
				
			}
			
			
			if ((OG1[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.13380E-03  1.28550E-07  ;  Og1-Hi+1\n", OG1[i]+1,  H[i+1]+1);
				
			}
			
			
			
			if (i>idx2res[m])
			{
				
				if ((OG1[i]>=0) && (O[i-1]>=0))
				{  
					npair++;fprintf(ftmp,"%d    %d  1  3.87420E-03   1.50095E-06 ;  Og1-Oi-1 \n", OG1[i]+1,  O[i-1]+1);
					
				}
				
				
				
			}
			
			
			
			if ((strcmp(resnm[i],"ASP")==0)||(strcmp(resnm[i],"ASN")==0)||(strcmp(resnm[i],"ASPH")==0))
			{
				
				if ((CG[i]>=0) && (H[i]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  0.0     8.5E-07  ;  CG-Hi asx\n", CG[i]+1,  H[i]+1);
					
				}
				if ((CG[i]>=0) && (H[i+1]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  0.0     8.5E-07  ;  CG-Hi+1 asx\n", CG[i]+1,  H[i+1]+1);
					
				}
				
				if ((CG[i]>=0) && (NH[i+1]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  1.91985E-03   2.06142E-06  ;  CG-Ni+1 asx \n", CG[i]+1,  NH[i+1]+1);
					
				}
				
				if ((ND2[i]>=0) && (NH[i]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  1.30345E-03   9.50217E-07  ;  ND2-Ni asx \n", ND2[i]+1,  NH[i]+1);
					
				}
				
				if ((ND2[i]>=0) && (NH[i+1]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1  2.76211E-03   4.26692E-06  ;  ND2-Ni+1 asx \n", ND2[i]+1,  NH[i+1]+1);
					
				}
				
				
				
				if (strcmp(resnm[i],"ASN")==0)
				{
					if ((HD21[i]>=0) && (O[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  1.275068E-03  2.139210E-08  ;  HD21-Oi asx\n", HD21[i]+1,  O[i]+1);
						
					}
					if ((OD1[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08  ;  OD1-Hi asx\n", OD1[i]+1,  H[i]+1);
						
					}
					
					
					if ((OD1[i]>=0) && (H[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08 ;  OD1-Hi+1 asx\n", OD1[i]+1,  H[i+1]+1);
						
					}
					
					
					
				}
				
				
				
				if ((strcmp(resnm[i],"ASP")==0)||(strcmp(resnm[i],"ASPH")==0))
				{
					
					
					if ((OD1[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08 ;  OD1-Hi asx\n", OD1[i]+1,  H[i]+1);
						
					}
					
					if ((OD1[i]>=0) && (H[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08 ;;  OD1-Hi+1 asx\n", OD1[i]+1,  H[i+1]+1);
						
					}
					
					if ((OD2[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08 ;;  OD2-Hi asx\n", OD2[i]+1,  H[i]+1);
						
					}
					
					
					if ((OD2[i]>=0) && (H[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.275068E-03  2.139210E-08 ;;  OD2-Hi+1 asx\n", OD2[i]+1,  H[i+1]+1);
						
					}
					
					
					if ((CG[i]>=0) && (O[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  7.41398E-03   2.29029E-06 ;;  CG-Oi asx\n", CG[i]+1,  O[i]+1);
						
					}
					
					
					
					
				}
				
				
				
				
				
			}
			
			else
				
				if ((strcmp(resnm[i],"PHE")==0)||(strcmp(resnm[i],"TYR")==0)||
					(strcmp(resnm[i],"TRP")==0)||(strcmp(resnm[i],"HIS")==0)||(strcmp(resnm[i],"HISE")==0)||(strcmp(resnm[i],"HISH")==0))
				{
					
					if ((CG[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  6.53184E-04  4.76171E-07  ;  CG-Hi aro\n", CG[i]+1,  H[i]+1);
						
					}
					if ((CG[i]>=0) && (H[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  6.53184E-04  4.76171E-07  ;  CG-Hi+1 aro\n", CG[i]+1,  H[i+1]+1);
						
					}
					
					
					if ((CG[i]>=0) && (NH[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1  1.91985E-03   2.06142E-06  ;  CG-Ni+1 aro \n", CG[i]+1,  NH[i+1]+1);
						
					}
					if ((CD1[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d 1  6.53184E-04  4.76171E-07  ;  Hi-CD1 aro\n", CD1[i]+1,  H[i]+1);
						
					}
					
					if ((ND1[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d 1  6.53184E-04  4.76171E-07  ;  Hi-ND1 aro\n", ND1[i]+1,  H[i]+1);
						
					}
					
					
					if ((CD2[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d 1  6.53184E-04  4.76171E-07  ;  Hi-CD2 aro\n", CD2[i]+1,  H[i]+1);
						
					}
					
					
					if ((CD1[i]>=0) && (NH[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.30345E-03   9.50217E-07  ;  CD1-Ni  aro\n", CD1[i]+1,  NH[i]+1);
						
					}
					
					if ((ND1[i]>=0) && (NH[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.30345E-03   9.50217E-07  ;  ND1-Ni  aro\n", ND1[i]+1,  NH[i]+1);
						
					}
					
					
					
					if ((CD2[i]>=0) && (NH[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  1.30345E-03   9.50217E-07  ;  CD2-Ni  aro\n", CD2[i]+1,  NH[i]+1);
						
					}
					
					
					
					if (i>idx2res[m])
					{
						
						if ((CD1[i]>=0) && (O[i-1]>=0))
						{
							npair++;fprintf(ftmp,"%d    %d   1     1.06354E-03  6.32621E-07  ;   CD1-Oi-1  aro\n", CD1[i]+1,  O[i-1]+1);
							
						}
						
						
						
					}
					
					if (i>idx2res[m])
					{
						
						if ((CD2[i]>=0) && (O[i-1]>=0))
						{
							npair++;fprintf(ftmp,"%d    %d    1     1.06354E-03  6.32621E-07  ;   CD2-Oi-1  aro\n", CD2[i]+1,  O[i-1]+1);
							
						}
						
						
						
					}
					
					if ((CD1[i]>=0) && (NH[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  2.76211E-03  4.26692E-06  ;  CD1-Ni+1 aro\n", CD1[i]+1,  NH[i+1]+1);
						
					}
					
					if ((ND1[i]>=0) && (NH[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  2.76211E-03  4.26692E-06  ;  ND1-Ni+1 aro\n", ND1[i]+1,  NH[i+1]+1);
						
					}
					
					
					
					if ((CD2[i]>=0) && (NH[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d   1  2.76211E-03  4.26692E-06  ;  CD2-Ni+1 aro\n", CD2[i]+1,  NH[i+1]+1);
						
					}
					
				}
				else
				{
					if ((CG[i]>=0) && (H[i]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1    1.74675E-03   1.70645E-06      ; H-CG H 0.2 0.24 nm rep\n", CG[i]+1,  H[i]+1);
						
					}
					if ((CG[i]>=0) && (H[i+1]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1    1.74675E-03   1.70645E-06      ; H-CG H 0.2 0.24 nm rep\n", CG[i]+1,  H[i+1]+1);
						
					}
					
					
					
					
					
				}
			
			
			if ((CG1[i]>=0) && (H[i]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.06354E-03  6.32621E-07      ; H-CG H 0.2 0.24 nm rep\n", CG1[i]+1,  H[i]+1);
				
			}
			if ((CG1[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.06354E-03  6.32621E-07      ; H-CG H 0.2 0.24 nm rep\n", CG1[i]+1,  H[i+1]+1);
				
			}
			
			
			
			
			if ((CG2[i]>=0) && (H[i]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.06354E-03  6.32621E-07      ; H-CG H 0.2 0.24 nm rep\n", CG2[i]+1,  H[i]+1);
				
			}
			if ((CG2[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1  1.06354E-03  6.32621E-07      ; H-CG H 0.2 0.24 nm rep\n", CG2[i]+1,  H[i+1]+1);
				
			}
			
			
			
			
			
			if ((NH[i]>=0) && (C[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1   3.37244E-03   5.68668E-06 ; 0.5 0.345 nm\n", NH[i]+1,  C[i+1]+1);
				
			}
			
			
			if ((O[i]>=0) && (O[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1   4.13270E-03   5.33725E-06 ; Oi-Oi+1 \n", O[i]+1,  O[i+1]+1);
				
			}
			
			
			
			
			
			
			if ((NH[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1   4.42E-04   1.64E-08  ;   Ni---Hi+1 \n", NH[i]+1,  H[i+1]+1);
				
			}
			
			
			
			
			
			if ((NH[i]>=0) && (O[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1   2.87870E-03  2.31737E-06 ; 0.894 0.305 nm\n", O[i+1]+1, NH[i]+1);
				
			}
			
			if ((CA[i]>=0) && (O[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1   2.26146E-03   3.19637E-06  ; 0.6  0.335 nm\n", CA[i]+1,  O[i+1]+1);
				
			}
			
			if ((C[i]>=0) && (O[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1    1.91985E-03   2.06142E-06  ; 0.6  0.32 nm\n", C[i]+1,  O[i+1]+1);
				
			}
			
			if ((O[i]>=0) && (C[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1    9.29809E-04   3.60227E-07  ; Oi-Ci+1\n", O[i]+1,  C[i+1]+1);
				
			}
			
			
			if ((CA[i]>=0) && (C[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1    2.49170E-03   6.92920E-06  ;  0.45  0.375 nm\n", CA[i]+1,  C[i+1]+1);
				
			}
			
			if ((CB[i]>=0) && (H[i+1]>=0))
			{
				npair++;fprintf(ftmp,"%d    %d  1    1.74675E-03   1.70645E-06      ; H-Cb H 0.2 0.24 nm rep\n", CB[i]+1,  H[i+1]+1);
				
			}
			
			
			
			if (i<nres-2)
			{
				
				if ((NH[i]>=0) && (C[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    3.37244E-03   5.68668E-06 ; 0.671 0.345 nm\n", NH[i]+1,  C[i+2]+1);
					
				}
				
				
				
				
				if ((NH[i]>=0) && (O[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.87870E-03  2.31737E-06 ; 0.894 0.305 nm\n", O[i+2]+1, NH[i]+1);
					
				}
				
				if ((CA[i]>=0) && (O[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.26146E-03   3.19637E-06   ; 0.6  0.335 nm\n", CA[i]+1,  O[i+2]+1);
					
				}
				
				if ((C[i]>=0) && (O[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    1.91985E-03   2.06142E-06   ; 0.6  0.32 nm\n", C[i]+1,  O[i+2]+1);
					
				}
				
				if ((C[i]>=0) && (NH[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1     3.37244E-03   5.68668E-06 ; 0.671 0.345 nm\n", C[i]+1,  NH[i+2]+1);
					
				}
				
				
				
				if ((O[i]>=0) && (CA[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.26146E-03   3.19637E-06  ; 0.6  0.335 nm\n", CA[i+2]+1, O[i]+1);
					
				}
				
				if ((O[i]>=0) && (NH[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    1.10468E-03   3.41254E-07 ; Oi-Ni+2\n", O[i]+1,  NH[i+2]+1);
					
				}
				
				
				if ((CA[i]>=0) && (CA[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.81500E-03   9.90525E-06 ; 0.45  0.39 nm\n", CA[i]+1,  CA[i+2]+1);
					
				}
				
				
				if ((CA[i]>=0) && (C[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.49170E-03   6.92920E-06  ;  0.45  0.375 nm\n", CA[i]+1,  C[i+2]+1);
					
				}
				
				
				if ((C[i]>=0) && (CA[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.49170E-03   6.92920E-06  ;  0.45  0.375 nm\n", C[i]+1,  CA[i+2]+1);
					
				}
				
				
				if ((C[i]>=0) && (C[i+2]>=0))
				{
					npair++;fprintf(ftmp,"%d    %d  1    2.17678E-03   4.73838E-06  ;  0.45  0.36 nm\n", C[i]+1,  C[i+2]+1);
					
				}
				
				
				
				if (i>idx2res[m])
				{
					if ((C[i-1]>=0) && (O[i+2]>=0))
					{
						npair++;fprintf(ftmp,"%d    %d  1     1.91985E-03   2.06142E-06   ; 0.6  0.32 nm\n", C[i-1]+1,  O[i+2]+1);
						
					}
				}
				
				
				
				
				
				
			}
			
			
			
			
			
			
			
			
			
			
			
		}
		
		
		
	}
	
	
	fclose (ftmp);

	//use info in temp.xvg
	int *pair1=malloc(npair*sizeof(int));	
	int *pair2=malloc(npair*sizeof(int));
	int *pairtyp=malloc(npair*sizeof(int));
	ftmp= fopen("temp.xvg","r");

	for (i=0;i<npair;i++)
	{
	 int a1, a2;
	 float v1, v2;
	 int isfound=0;
	 get_line (ftmp,buffer);
	 sscanf (buffer,"%d %d %*d %f %f\n", 
		&a1, &a2, &v1, &v2);
	pair1[i]=a1;
	pair2[i]=a2;
	
	  for (j=0;j<nparam;j++)
	   if (iseq(v1,vdwParam[j][0],1e-4)&&iseq(v2,vdwParam[j][1],1e-4))
	     {isfound=1;pairtyp[i]=j+1;break;}
	  if (!isfound)
	  {
		if (nparam==MAXINTTYP) {printf("# of modified pair types is larger than the limit %d.",MAXINTTYP);exit(1);}
		vdwParam[nparam][0]=v1;
		vdwParam[nparam][1]=v2;
		pairtyp[i]=nparam+1;
		nparam++;
          }
	}	

	fclose(ftmp);
//add here the important information about lipid exclusion
//start with islipid[i] 0 1 for all atoms
//nlipatom
       fclose(fp);
if (nlipatom)
{
 #define MAXCON 10
   int **lipcon = malloc(nlipatom*sizeof(int*));
   int *lipconcnt = malloc(nlipatom*sizeof(int));
   int *lip2aidx = malloc(nlipatom*sizeof(int));
   for(i=0;i<psfnatom;i++)
    if (islip[i]>=0)
      lip2aidx[islip[i]]=i+1;
   for (i=0;i<nlipatom;i++)
    {
      lipcon[i]=malloc(MAXCON*sizeof(int));
       lipconcnt[i]= 0 ;
//       for (j=0;j<MAXCON;lipcon[i][j++]=-1);
   }

   fp=fopen(argv[1],"r");
//scan bond info from psf file
  int found=0;
   while (!feof(fp))
  {
     get_line(fp,buffer);
     trcspace(buffer);
     if (parcnt(buffer)>1)
     {
      parstr(buffer,2,buf);
       if (!strncmp(buf,"!NBO",4))
         {
           found=1; break;
        }
      }    
  };

  if (!found)
  {
    printf("Error: No bond connections are found in the PSF.\n");
    exit(1);
  }

   while (!feof(fp))
  {
     get_line(fp,buffer);
     trcspace(buffer);
     if (parcnt(buffer)>1)
     {
      parstr(buffer,2,buf);
       if (!strncmp(buf,"!N",2))
         {
           found=1; break;
        }
    //load bond idx pair 
     int nidx=parcnt(buffer);
       for (i=1;i<=nidx;i+=2)
        {
           int idxa, idxb;
           parstr(buffer,i,buf);
           sscanf(buf,"%d",&idxa);
           parstr(buffer,i+1,buf);
           sscanf(buf,"%d",&idxb);

          if (islip[idxa-1]>=0&&islip[idxb-1]>=0)
          {
             int i_a  =  islip[idxa-1];
             int i_b  =  islip[idxb-1];
             if (lipconcnt[i_a]>=MAXCON||lipconcnt[i_b]>=MAXCON)
             {
               printf("Error: Too many connections to atom \n");
               exit(1);
            }
             lipcon[i_a][(lipconcnt[i_a])++]=idxb;
             lipcon[i_b][(lipconcnt[i_b])++]=idxa;
          }


        }

      }

    
  };   


//search for 13 interaction
char *mark= malloc(nlipatom*sizeof(char));
int *queue= malloc(nlipatom*sizeof(int));
int *gene = malloc(nlipatom*sizeof(int));

 struct LLpair {
 int ia;
 int ib;
 struct LLpair * next; 
}; 

struct LLpair ts, *tn;
 tn = &ts;
 tn->next=NULL;
int nnewpair=0;

for (i=0;i<nlipatom;i++)
{
  for (j=0;j<nlipatom;mark[j++]=0);
  int head=0, reer=0;
  queue[head]=i;
  mark[i]=1;
  gene[head]=1;
  while (head<=reer)
  {
     int _gene = gene[head];
     int _ci = queue[head];
     if (_gene==3)
     {
       if (lip2aidx[i]<lip2aidx[_ci])
       {
         //add pair
           tn->next = malloc(sizeof(struct LLpair));
           tn= tn->next;
           tn->ia =   lip2aidx[i];
             tn->ib =   lip2aidx[_ci];
            tn->next=NULL;
             nnewpair++;
         }
       head++; continue;
     }
 
     for (j=0;j<lipconcnt[_ci];j++)
     {
       int _lid = islip[lipcon[_ci][j]-1];
       if (mark[_lid]) continue;
       mark[_lid] =1;
       reer++;
       queue[reer] = _lid;
       gene[reer] = _gene +1;       
       
     }
     head++;
  }  
}
//all new pair stored in ts (LLpair)

 
 int *t_p1=malloc((npair+nnewpair)*sizeof(int));
 int *t_p2=malloc((npair+nnewpair)*sizeof(int));
 for (i=0;i<npair;i++)
 {
   t_p1[i]=pair1[i];
   t_p2[i]=pair2[i];
 }

 tn = ts.next;
 while (tn!=NULL)
 {
   struct LLpair *_t = tn;
   t_p1[npair]=tn->ia;
   t_p2[npair]=tn->ib;
   //printf("%d %d\n",tn->ia, tn->ib);
   npair++;
   tn = tn->next;
   free(_t);
 }

 free(pair1);
 free(pair2);
 pair1 = t_p1;
 pair2 = t_p2;
 t_p1=t_p2=NULL;

   fclose(fp);

   for(i=0;i<nlipatom;free(lipcon[i++]));
   free(lipcon);
   free(lipconcnt);
   free(lip2aidx);
   free(mark);
   free(queue);
   free(gene);
}



#ifdef DEBUG
  psfnatom = 5;
  npair= 3;
  pair1[0] = 4; pair2[0]=3;
  pair1[1] = 5; pair2[1]=3;
  pair1[2] = 3; pair2[2]=2;

#endif
//sort pair1[i] pair2[i]

 for (i=0;i<npair;i++)
  if (pair1[i]>pair2[i])
  {
    int ii = pair1[i];
    pair1[i]=pair2[i];
    pair2[i]= ii;
   }


  for (i=0;i<npair-1;i++)
   for (j=i+1;j<npair;j++)
    if (pair1[i]>pair1[j])
    {
    int ii = pair1[i];
    pair1[i]=pair1[j];
    pair1[j]= ii;


        ii = pair2[i];
    pair2[i]=pair2[j];
    pair2[j]= ii;


    }
   else if ((pair1[i]==pair1[j])&&
   (pair2[i]>pair2[j]))
   {
   int  ii = pair2[i];
    pair2[i]=pair2[j];
    pair2[j]= ii;

    }

  //now for all pairs like
  // 1  1   1   2  2  3  3  3  4   4  4 ..
  // 2  4   5   3  6  5  6  7  5   8  9 .. 


#if 0
	ftmp=fopen(argv[2],"r");
	fp=fopen(argv[4],"w");
	while (!feof(ftmp))
	{
		get_line(ftmp,buffer);
		fprintf(fp,"%s\n",buffer);	
        }
	fclose(ftmp);

	fprintf(fp,"\n");
	fprintf(fp,"MODP\n");
	for (i=0;i<nparam;i++)
	fprintf (fp,"%.6g  %.6g  %.6g  %.6g  ! %d\n", 
		vdwParam[i][1]/4.184*1.e+12,vdwParam[i][0]/4.184*1.e+6,
 		vdwParam[i][1]/4.184*1.e+12,vdwParam[i][0]/4.184*1.e+6,i+1);
	
	fprintf(fp,"\n");
	fclose(fp);
#endif
	ftmp=fopen(argv[1],"r");
	fp=fopen(argv[2],"w");
int showatom=0;
  int wrapcnt=0;	
int turnoff=0;
int lastres, suf_id;
	while (!feof(ftmp))
	{
		get_line(ftmp,buffer);
		strcpy(buf,buffer);
		trcspace(buf);
		if (parcnt(buf)>1)
		{
		  char ct[512];
		  parstr(buf,2,ct);

                if (!strcmp(ct,"!NATOM")) 
                 { 
                    showatom=1;turnoff=1; 
                     lastres=-1; suf_id=-1;
		fprintf(fp,"%s\n",buffer);
                    continue;}

		  if (!strcmp(ct,"!NNB"))
		  {
			fprintf(fp,"%8d !NNB\n",npair);
                  //output pairs
//first part of NNB
                   for (i=0;i<npair;i++)
                    {
                     fprintf(fp,"%8d",pair2[i]);
                     wrapcnt++;
                     if (wrapcnt==8) 
                     {
                        wrapcnt=0; 
                        fprintf(fp,"\n");
                      }
                    }
                    if (wrapcnt) fprintf(fp,"\n");
fprintf(fp,"\n");
//second part of NNB
int idx_pair1 =0;
wrapcnt=0;
 for (i=0;i<psfnatom;i++)
{
  int c_idx = i+1;
  if (idx_pair1<npair)
   if (c_idx == pair1[idx_pair1])
   {
      for (j=idx_pair1;j<npair;j++)
       if (pair1[j]>c_idx) break;
      idx_pair1 = j;
   }
  fprintf(fp,"%8d",idx_pair1);
                     wrapcnt++;
                     if (wrapcnt==8)
                     {
                        wrapcnt=0;
                        fprintf(fp,"\n");
                      }

}

if (wrapcnt) fprintf(fp,"\n");

fprintf(fp,"\n");

              turnoff=1;
                  continue;

//			for (i=0;i<npair;i++)
//			fprintf (fp,"%8d%8d%8d\n",pair1[i],pair2[i],pairtyp[i]);
//			fprintf(fp,"\n");
		  }
		}
         if (turnoff)
         {
           //fprintf(fp,"off %d\n", strlen(buffer)); 
           for (i=0;i<(int)(strlen(buffer))-1;i++)
             if (buffer[i]=='!'&&buffer[i+1]=='N')
              {
            //       fprintf(fp,"on %d\n",i);
                  showatom=0;
                    turnoff=0; break;
               }
         }

        if (showatom&&parcnt(buf)>1)
        {
           char ct[512];
            int resid;
             
             parstr(buf,3,ct);
            sscanf(ct,"%d", &resid);
            
            parstr(buf, 6, ct);
            int lc=(int)strlen(ct);

            int st= strindex(buffer, ct, 0);
            if (resid!=lastres)
            {
                suf_id++;
                if (suf_id==2*NFORWARD+1) suf_id=0;
               lastres=resid;
            }
           char resnm[512];
           parstr(buf, 4, resnm);
           int ispro=0;
           for (i=0;i<NAMINO;i++)
            if (!strcmp(resnm, aanm[i]))
              {ispro=1;break;}
          if (ispro)
         {
            ct[lc]='x';ct[lc+1]='A'+(char)suf_id; ct[lc+2]=0;            
          }
  
            buffer[st] = 0;
            char *suf = buffer + st + (int)strlen(ct);
            fprintf(fp,"%s%s %s\n", buffer,ct,suf);
        }
        else if (showatom) fprintf(fp,"\n");
        
         if (!turnoff)	
		fprintf(fp,"%s\n",buffer);	
	}

	fclose(ftmp);
	fclose(fp);



	free (pair1); free(pair2); free(pairtyp);

	//
	free(C);
	free(CA);
	free(NH);
	free(CB);
	free(O);
	free(H);
	free(CG1);
	free(CG2);
	free(CG);
	free (OG);
	free(OE);
	free(NE);
	free(AUX);
	free(CD1);
	free(CD2);
	free (OG1);
	free(OD1);
	free(OD2);
	free(HD1);
	free(ND1);
	free(ND2);
	free(HD21);
	free(HD22);
	free (idx2res);
	for (i=0;i<nproatom+1;i++) free(resnm[i]);
	free (resnm);
	
	
	
	
	for (i=0;i<psfnatom;i++)
	{
		free(psfsegnm[i]);
		free(psfresnm[i]);
		free(psfatomnm[i]);
		free(psfatomtype[i]);
        free(psfnewnm[i]);
		
	}
	
	free (psfnewnm);
	free(psfsegnm);
	free(psfresnm);
	free(psfatomnm);
	free(psfatomtype);
	free(psfresid);
        free(islip);	
	return 0;
	
}