autotors.c

c++编写的并行拉马克遗传算法的程序。实现分析对接程序

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>
#include <sys/times.h>
#include <ctype.h>
#include <malloc.h>
#define  MAX_BNDFLD  20  /*length of bnd listing*/
#define  MAX_CONNECTS  6
#define  MAX_CYCLES  10
#define  MAX_MEMBERS  20 /*largest possible cycle */
#define  MAX_RECORD  200
#define  LINE_LEN  120
#define  MAX_TORS  64
#define  equal(a,b,n) (strncmp(a,b,(n)) == 0)
#define  UNKNOWNREC  0
#define  ATOMREC  1
#define  BONDREC  2
#define  ANCHORREC  3
#define  RBONDREC  4
#define  MOLECULEREC  5

/*----------------------------------------------------------------------------*/
typedef struct ATOM_NODE
{
	int id; /*atom number is its id*/
	int new_id; /*re-numbered atom  id*/
	int root; /*flag to mark root-atoms*/
	struct ATOM_NODE *prev[MAX_CONNECTS]; /*parent*/
	struct ATOM_NODE *next[MAX_CONNECTS]; /*children*/
	int prev_num; /*number of previous links*/
	int next_num; /*number of next links*/
} ATOM_NODE;

ATOM_NODE * head;

typedef struct TORS_ANGLE
{
	int id1; /*one of the two old ids involved*/
	int id2; /*the other one*/
	int num; /*counter for the tors*/
	char status;/*'I' for inactive torsions; 'A' for active ones*/
	char atom1_rec[6],atom2_rec[6]; /*atom name: element+number */
	struct TORS_ANGLE *next; /*the next in the linked list of torsions*/
} TORS_ANGLE;

TORS_ANGLE *  head_tors, *next_tors, * temphead_tors;


int debug;
int *find_key_list,find_key_num;
int hflag; /*optional merging on non-polar hydrogens with their carbons*/
int mflag; /*used for mol2 format bnd input*/
int aflag; /*when present, amide bond torsions are NOT allowed*/
int cflag; /*causes program to output number of connections for each atom*/
int bflag; /*disallows torsions in peptide backbone*/
int Mflag; /*with mflag causes automatic root selection via ROTATABLE_BOND+ANCHOR sections */
int rflag; /*enable automatic selection of atom closest to the center as root*/
int Aflag; /*program tests if rings aromatic and replaces aromatic C with A in output*/
int offset;/*set to 70 or 55 depending on which pdbq format input*/
int nbranch;/*number of current open branches;used as flag for root expansion*/
int ntors;/*number of tors:used to check whether MAX_TORS exceeded*/
int stors;/*number of selected tors*/
int torsdof;/*number of tors used for force field calculations*/
int rootnum; /*the number of atoms defined to be in the root*/
int rootlist[MAX_RECORD];
int cycles[MAX_CYCLES][MAX_MEMBERS];/*to store cycles detected*/
int cycle_size[MAX_CYCLES];/*number of members of each cycle*/
int cycle_count; /*count number of cycles detected*/
int traverse_num; /*flag to mark which time to traverse tree structure*/
int nonpolarH; /*counter for non-polar H's merged with their carbon atoms*/
char *atmptr;/*data structure pointer for pdbq records*/
int natom,outatom;/*number of atoms found in input file;number of atoms output*/
FILE *outfile, *file_ptr;
void traverse_cycle();


char atmtype(charptr, linenum)
char * charptr;
int linenum;

{
	char c12,c13;

	c12= *(charptr+LINE_LEN*linenum+12);
	c13= *(charptr+LINE_LEN*linenum+13);
/*fprintf(stderr, "c12=%c, c13=%c\n",c12,c13);*/
	if(c12==' ') {
		if(c13=='C') {
			if (*(charptr+LINE_LEN*linenum+14)=='A'||
			*(charptr+LINE_LEN*linenum+14)=='a'){
				return('a');
/*ASSUME CA not involved in amide bond ever so don't need to return plain "c"*/
			};
			if (*(charptr+LINE_LEN*linenum+14)=='B'||
			*(charptr+LINE_LEN*linenum+14)=='b'){
				return('b');
/*detect CB not involved in amide bond ever so don't need to return plain "c"*/
			};
			if (*(charptr+LINE_LEN*linenum+14)=='L'||
			*(charptr+LINE_LEN*linenum+14)=='l'){
				return('X');
			}else return(c13);
		}else return(c13);
	}
	return(c12);	
}

ATOM_NODE * get_new_node(iden)
int iden;
{
	int i;
	ATOM_NODE * temp;
	temp=(ATOM_NODE *) malloc(sizeof(ATOM_NODE));
	for(i=0;i<MAX_CONNECTS;i++){
		temp->prev[i]=NULL;
		temp->next[i]=NULL;
	}
	temp->id=0;
	temp->new_id=0;
	temp->root=0; /*assume not a root until told otherwise*/
	temp->next_num=0;
	temp->prev_num=0;
	temp->id=iden;
	return(temp);
}


ATOM_NODE * first(key1,key2)
int key1;
int key2;

{    
	ATOM_NODE  *current;
	head=get_new_node(key1);
	current=get_new_node(key2);
	head->next[0]=current;
	head->next_num++;
	current->prev[0]=head;
	current->prev_num++;
	return(head);
}

int  on_list(key,list,index)
int key, * list, index;
{
	int i,ans;
	ans=0;
	for(i=0;i<index;i++){
		if(*(list+i)==key) {
			ans=1;
		};
	};
	return(ans);
}

ATOM_NODE * check_node(key,node,clist,cindex)
int key;
ATOM_NODE * node;
int * clist,* cindex;

{
	int i;
	ATOM_NODE * ans;

	ans=NULL;
	*(clist+*cindex)=node->id;
	*cindex=*(cindex)+1;

	/*just have added current node id to visited list and incremented counter*/
	if(node->id==key){
		return(node);
	};
	if(node->prev_num){
		for(i=0;i<node->prev_num;i++){
			if(!on_list((node->prev[i])->id,clist,*cindex)){
				ans=check_node(key,node->prev[i],clist,cindex);
				if(ans) return(ans);
			};
		};
	};
	if(node->next_num){
		for(i=0;i<node->next_num;i++){
			if(!on_list((node->next[i])->id,clist,*cindex)){
				ans=check_node(key,node->next[i],clist,cindex);
				if(ans) return(ans);
			};
		};
	};
	return(ans);
}

ATOM_NODE * find_key(key,node_num,node)
int key,node_num;
ATOM_NODE * node;

{
	int i;
	ATOM_NODE * ans;

        for(i=0;i<node_num;i++) *(find_key_list+i)=0;
        find_key_num=0;
	ans=NULL;
	ans=check_node(key,node,find_key_list, &find_key_num);
	return(ans);
}



ATOM_NODE * mark_node(key,node,mlist,mindex,vlist,vindex)
int key;
ATOM_NODE * node;
int * mlist,*mindex,*vlist,*vindex;

{    
	int i;
	ATOM_NODE *ans;

	ans = NULL;
	/*mlist is list of atoms to be marked as root atoms*/
	*(mlist+*mindex)=node->id;
	*mindex=*mindex + 1;

	/*vlist is list of atoms visited, not necessarily to be marked*/
	*(vlist+*vindex)=node->id;
	*vindex=*vindex +1;

	if(node->id==key){
		return(node);
	};
	if(node->prev_num){
		for(i=0;i<node->prev_num;i++){
			if(!on_list((node->prev[i])->id,vlist,*vindex)){
				ans=mark_node(key,node->prev[i],mlist,mindex,
                                              vlist,vindex);
				if(ans) return(ans);
			};
		};
	};
	if(node->next_num){
		for(i=0;i<node->next_num;i++){
			if(!on_list((node->next[i])->id,vlist,*vindex)){
				ans=mark_node(key,node->next[i],mlist,mindex,
                                              vlist,vindex);
				if(ans) return(ans);
			};
		};
	};
	if(!ans){
		*mindex=*mindex-1;
		*(mlist+*mindex)=0;
	};
	return(ans);
}

/*When a cycle is detected, mark all atoms in it as ROOT atoms*/
ATOM_NODE * mark_atom(key,con_num,node)
int key,con_num;
ATOM_NODE * node;

{    
	int i,j,k,v_num,m_num;
	int p,cindex,cadded;
	int * visited,*marked;
	ATOM_NODE * ans, *temp;

	/*first get space for the list of atoms to be visited at most con_num*/
	visited = (int*) calloc (con_num,sizeof(v_num));
	marked = (int*) calloc (con_num,sizeof(v_num));
	ans=NULL;
	v_num=0;
	m_num=0;
	cadded=0;

	ans=mark_node(key,node,marked, &m_num,visited,&v_num);

	for(i=0;i<m_num;i++){
	   temp=find_key(*(marked+i),con_num,head);
	   if(temp&&temp->root==1){
	      for(j=1;j<=cycle_count;j++){
	         for(k=0;k<cycle_size[j];k++){
	            if(cycles[j][k]==temp->id){
		       cycle_count--;
		       cindex=cycle_size[j];
		       for(p=0;p<m_num;p++){
			  if(!on_list(*(marked+p),cycles[j],cindex+cadded)){
			     cycles[j][cindex+cadded]=*(marked+p);
			     temp=find_key(*(marked+p),con_num,head);
			     temp->root=1;
			     cadded++;
			  };
                /*need to add to existing cycle+MARK ALL roots*/
		/*need to update cycle_size[cycle_count]*/
		       };
		       cycle_size[j]=cycle_size[j]+cadded;
/*having added the new members to an existing cycle + updated count,done*/
		       return;
		    };
		 };
	      };
	   printf("error in root designation of %d atom\n",temp->id);
	   exit(1);
	   };
	};
	/*just dealt with adding to pre-exixting cycle*/
	for(i=0;i<m_num;i++){
	   temp=find_key(*(marked+i),con_num,head);
	   if(temp){
	      temp->root=1;
	      cycles[cycle_count][i]=temp->id;
	   };
	   if(!temp)printf("error: marking roots-couldn't find %d\n",temp->id);
        };
/*at this point have added a new cycle to list and must record its size*/
	cycle_size[cycle_count]=m_num;
	free(visited);
	free(marked);
	return(ans);
}

int* is_cycle(node)
ATOM_NODE * node;

{	int i,j,k,n,m;
	int *newlist,newindex,entry;

for(j=1;j<=cycle_count;j++){
	for(k=0;k<cycle_size[j];k++){
		if(cycles[j][k]==node->id){
			if(node->root){
       			newlist= (int*) calloc(cycle_size[j],sizeof(newindex));
			newindex=0;
   			for(i=cycle_size[j];i<2*cycle_size[j];i++){
			*(newlist+newindex)=cycles[j][(k+i)%cycle_size[j]];
			newindex++;
 			};/*for all cycle entries*/
                        return(newlist);
                       	} else{
                      	printf("ERR:root designation of %d\n",node->id);
                       	exit(1);
			};/*if-else*/
		};
	};/*for k's*/
};/*for j's*/
return (NULL);
}

int num_cycle(node)
ATOM_NODE * node;

{ 	int j,k;
for(j=1;j<=cycle_count;j++){
	for(k=0;k<cycle_size[j];k++){
		if(cycles[j][k]==node->id){ 
			return(j);
		};
	};
 };
return(0);
}

int eq_cycle(key1,key2)
int  key1,key2;

{
        int j,k;
        int num1,num2;
        num1=-1;
        num2=-1;
        for(j=1;j<=cycle_count;j++){
                for(k=0;k<cycle_size[j];k++){
                        if(cycles[j][k]==key1){
                                num1=j;
                        };
                        if(cycles[j][k]==key2){
                                num2=j;
                        };
                };
         };
        if((num1<0)||(num2<0)){
                 return(0);
        };
        if(num1!=num2){
                 return(0);
        }
        return(1);
}


TORS_ANGLE * get_new_tors(id1,id2)
int id1;
int id2;

{    
	int j;
	char *tptr1,*tptr2;
	TORS_ANGLE * temp;
	temp=(TORS_ANGLE *) malloc (sizeof(TORS_ANGLE));

	if(!id1||!id2){
		printf("error in get_new tors\n");
		exit(1);
	};
	temp->id1=id1;
	temp->id2=id2;
	temp->num=ntors;
	temp->status='A';

	for(j=0;j<6;j++){
		temp->atom1_rec[j]=0;
		temp->atom2_rec[j]=0;
	};
	j=0;

	if(*(atmptr+LINE_LEN*id1+13+j)==' '){
		printf("error in get_new tors on id1\n");
		exit(1);
	};
/*copy space 12 also*/
/*	for(j=0;j<4;j++)
		temp->atom1_rec[j]=*(atmptr+LINE_LEN*id1+13+j);*/
		
	for(j=0;j<5;j++)
		temp->atom1_rec[j]=*(atmptr+LINE_LEN*id1+12+j);
	temp->atom1_rec[j]='\0';

	j=0;

	if(*(atmptr+LINE_LEN*id2+13+j)==' '){
		printf("error in get_new tors on id2\n");
		exit(1);
	};

	for(j=0;j<5;j++)
		temp->atom2_rec[j]=*(atmptr+LINE_LEN*id2+12+j);
	temp->atom2_rec[j]='\0';

	tptr1=temp->atom1_rec;
	tptr2=temp->atom2_rec;

	temp->next=NULL;
	return(temp);
}

int  check_tors(key1,key2)