autotors.c

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

	for(i=0;i<root_num;i++){
		temp_root=find_key(rootlist[i],con_num,head);
		if(!temp_root){
			printf("error in finding %d\n",rootlist[i]);
			exit(1);
		};

		if((!temp_root)||(!temp_root->root)) {
			fprintf(outfile, "%d atom not a root\n",temp_root->id);
			exit(1);
		};
		connect_num=temp_root->prev_num + temp_root->next_num;
		p_num=temp_root->prev_num;
		n_num=temp_root->next_num;
		switch(connect_num){
		case 0 : break;
		case 1 : 
if(temp_root->prev_num&&!on_list(temp_root->prev[0]->id,travlist,trav_index)){
    if(!eq_cycle(temp_root->prev[0]->id,temp_root->id)){
	if(traverse_num==1){
	     next_prev(temp_root->prev[0],travlist,&trav_index,nplist,np_index);
	     for(k=0;k<*np_index;k++){
	     *(nplist+k)=0;
	     };
	     *np_index=0;
	};
	visit_prev(temp_root->prev[0],travlist,&trav_index,nplist,np_index);
	break;
  };
};
if(temp_root->next_num&&!on_list(temp_root->next[0]->id,travlist,trav_index)){
    if(!eq_cycle(temp_root->next[0]->id,temp_root->id)){
	if(traverse_num==1){
	     prev_next(temp_root->next[0],travlist,&trav_index,nplist,np_index);
	     for(k=0;k<*np_index;k++){
	     *(nplist+k)=0;
	     };
	     *np_index=0;
	};
	visit_next(temp_root->next[0],travlist,&trav_index,nplist,np_index);
	break;
	};
};
break;

	case 2 : for(j=0;j<p_num;j++){
   	prev_con=temp_root->prev[j]->prev_num+temp_root->prev[j]->next_num;
   	if(!on_list(temp_root->prev[j]->id,rootlist,rootnum)&&!on_list(temp_root->prev[j]->id,travlist,trav_index)){
    if(!eq_cycle(temp_root->prev[j]->id,temp_root->id)){
	 if(traverse_num==1){
	     next_prev(temp_root->prev[j],travlist,&trav_index,nplist,np_index);
	     for(k=0;k<*np_index;k++){
	     *(nplist+k)=0;
	     };
	     *np_index=0;
	 };
	};
	if((traverse_num==1)&&(prev_con>1)){
	if(aflag && is_amide(temp_root)<0)printf("8:found amide prev %d\n",temp_root->id);
	else	
	if(!(bflag && check_back(temp_root,j,'p')))
	head_tors=add_tors(temp_root->id,temp_root->prev[j]->id,head_tors);
	};
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->prev[j]->id))
	nbranch++;
	if(traverse_num==2&&!nbranch&&!check_tors(temp_root->id,temp_root->prev[j]->id)){
		 if(!on_list(temp_root->prev[j]->id,rootlist,rootnum)&&!eq_cycle(temp_root->prev[j]->id,temp_root->id)){

		temp_root->prev[j]->root=1;
		rootlist[rootnum++]=temp_root->prev[j]->id;
		};
	};
        if(traverse_num==4&&check_tors(temp_root->id,temp_root->prev[j]->id))
fprintf(outfile,"BRANCH %d %d\n",temp_root->new_id,temp_root->prev[j]->new_id);
	visit_prev(temp_root->prev[j],travlist,&trav_index,nplist,np_index);
      	if(traverse_num==2&&check_tors(temp_root->id,temp_root->prev[j]->id))
	nbranch--;
      	if(traverse_num==4&&check_tors(temp_root->id,temp_root->prev[j]->id))
fprintf(outfile,"ENDBRANCH %d %d\n",temp_root->new_id,temp_root->prev[j]->new_id);
	};
 	};/*for all prevs*/
	for(j=0;j<n_num;j++){
	next_con=temp_root->next[j]->prev_num+temp_root->next[j]->next_num;
   	if(!on_list(temp_root->next[j]->id,rootlist,rootnum)&&!on_list(temp_root->next[j]->id,travlist,trav_index)){
        if(!eq_cycle(temp_root->next[j]->id,temp_root->id)){

	 if(traverse_num==1){
	     prev_next(temp_root->next[j],travlist,&trav_index,nplist,np_index);
	     for(k=0;k<*np_index;k++){
	     *(nplist+k)=0;
	     };
	     *np_index=0;
	 };
	};
	if(traverse_num==1&&(next_con>1)){
	if(aflag && is_amide(temp_root)>0)printf("9:found amide next %d\n",temp_root->id);
	else
	     if(!(bflag && check_back(temp_root,j,'n')))
	     head_tors=add_tors(temp_root->id,temp_root->next[j]->id,head_tors);
	};
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->next[j]->id))
	nbranch++;
	if(traverse_num==2&&!nbranch&&!check_tors(temp_root->id,temp_root->next[j]->id)){
if(!on_list(temp_root->next[j]->id,rootlist,rootnum)&&!eq_cycle(temp_root->next[j]->id,temp_root->id)){
		temp_root->next[j]->root=1;
		rootlist[rootnum++]=temp_root->next[j]->id;
		};
	};
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->next[j]->id))
fprintf(outfile,"BRANCH %d %d\n",temp_root->new_id,temp_root->next[j]->new_id);
	visit_next(temp_root->next[j],travlist,&trav_index,nplist,np_index);
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->next[j]->id))
	nbranch--;
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->next[j]->id))
fprintf(outfile,"ENDBRANCH %d %d\n",temp_root->new_id,temp_root->next[j]->new_id);
	};
  	};/*for all nexts*/
	break;

	default:for(j=0;j<p_num;j++){
   	prev_con=temp_root->prev[j]->prev_num+temp_root->prev[j]->next_num;
   	if(!on_list(temp_root->prev[j]->id,rootlist,rootnum)&&!on_list(temp_root->prev[j]->id,travlist,trav_index)){
        if(!eq_cycle(temp_root->prev[j]->id,temp_root->id)){
	 if(traverse_num==1){
	     next_prev(temp_root->prev[j],travlist,&trav_index,nplist,np_index);
     	     for(k=0;k<*np_index;k++){
		*(nplist+k)=0;
	     };
	     *np_index=0;
	 };
	};
	if((traverse_num==1)&&(prev_con>1)){
		if(aflag && is_amide(temp_root)<0)printf("10:found amide prev %d\n",temp_root->id);
		else
	        if(!(bflag && check_back(temp_root,j,'p')))
		head_tors=add_tors(temp_root->id,temp_root->prev[j]->id,head_tors);
	};
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->prev[j]->id))
	nbranch++;
	if(traverse_num==2&&!nbranch&&!check_tors(temp_root->id,temp_root->prev[j]->id)){
                if(!on_list(temp_root->prev[j]->id,rootlist,rootnum)&&!eq_cycle(temp_root->prev[j]->id,temp_root->id)){
		temp_root->prev[j]->root=1;
		rootlist[rootnum++]=temp_root->prev[j]->id;
		};
	};
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->prev[j]->id))
fprintf(outfile,"BRANCH %d %d\n",temp_root->new_id,temp_root->prev[j]->new_id);
	/*in any case visit the node*/
	visit_prev(temp_root->prev[j],travlist,&trav_index,nplist,np_index);
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->prev[j]->id))
	nbranch--;
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->prev[j]->id))
	fprintf(outfile,"ENDBRANCH %d %d\n",temp_root->new_id,temp_root->prev[j]->new_id);
	};
	}; /*for-prev*/
	for(j=0;j<n_num;j++) {
	next_con=temp_root->next[j]->prev_num+temp_root->next[j]->next_num;
   	if(!on_list(temp_root->next[j]->id,rootlist,rootnum)&&!on_list(temp_root->next[j]->id,travlist,trav_index)){
        if(!eq_cycle(temp_root->next[j]->id,temp_root->id)){
	 if(traverse_num==1){
	     prev_next(temp_root->next[j],travlist,&trav_index,nplist,np_index);
	     for(k=0;k<*np_index;k++){
		*(nplist+k)=0;
	     };
	     *np_index=0;
	 };
	};
	if(traverse_num==1&&(next_con>1)){
		if(aflag && is_amide(temp_root)>0)printf("11:found amide next %d\n",temp_root->id);
		else
	        if(!(bflag && check_back(temp_root,j,'n')))
   	        head_tors=add_tors(temp_root->id,temp_root->next[j]->id,head_tors);
	};
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->next[j]->id))
	nbranch++;
	if(traverse_num==2&&!nbranch&&!check_tors(temp_root->id,temp_root->next[j]->id)){
                if(!on_list(temp_root->next[j]->id,rootlist,rootnum)&&!eq_cycle(temp_root->next[j]->id,temp_root->id)){
		temp_root->next[j]->root=1;
		rootlist[rootnum++]=temp_root->next[j]->id;
		};
	};
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->next[j]->id))
fprintf(outfile,"BRANCH %d %d\n",temp_root->new_id,temp_root->next[j]->new_id);
	/*in any case visit the node*/
	visit_next(temp_root->next[j],travlist,&trav_index,nplist,np_index);
	if(traverse_num==2&&check_tors(temp_root->id,temp_root->next[j]->id))
	nbranch--;
	if(traverse_num==4&&check_tors(temp_root->id,temp_root->next[j]->id))
	fprintf(outfile,"ENDBRANCH %d %d\n",temp_root->new_id,temp_root->next[j]->new_id);
	};
	}; /*for-next*/
	break;
}; /*end of switch  on con_num statement*/
};/*for each rootlist entry*/
	free(nplist);
}

void main(argc,argv)

int argc;
char **argv;
/*----------------------------------------------------------------------------*/

{
	int m,k,j,i,ir,ii;
	int nrecord;
	int newroot;
	int *visited;
	int oldmin,newmin,oldmax;
	int sort_root[MAX_RECORD];
	int cycleans;
	int argindex;
	int mindex;
	int foundtors;
	char rootans,rootkey[6];
	char newans[3];
	char line[LINE_LEN], rec4[4];
	int atom_1,atom_2,con_num;   /*number of connections input*/
	int treenum; /*number of temporary trees in forest*/
	int treecount; /*number of temporary trees in forest*/
	int treetry; /*number of times tried to attach temporary trees*/
	char record[MAX_RECORD][LINE_LEN];
	char pdbq_record[MAX_RECORD][LINE_LEN];
	ATOM_NODE   *temphead[MAX_RECORD], *temp1, *temp2;
	TORS_ANGLE * temp, *t2;
	int anum;
	int resnum;
	char pdbaname[7];
	char resname[80];
	float x,y,z;
	float charge;
	float xcen;
	float ycen;
	float zcen;
	float newx;
	float newy;
	float newz;
	float diffx;
	float diffy;
	float diffz;
	int offsetx;
	int offsety;
	int offsetz;
	int nearest;
	int ij;
	float neardist;
	float newdist;

	FILE *connect_file, *pdbq_file;

	typedef struct BOND_STRUCT{
		int atom1;
		int atom2;
	} BOND;
	int natoms,nanchors,nrbonds,bnum,stat,ring,recordtype,linenum;
	BOND *bp, *bndptr;
	int rbond[MAX_TORS];
	typedef struct EQV_STRUCT{
		char sybtyp[6];
		char autotyp[2];
	} EQVREC;
	EQVREC *tp,eqvtable[40];
	int eflag,neqvs;
	FILE *eqv_file;
	char atyp[6];

	debug=0;
	ntors=0;
	nbranch=0;
	con_num=0;
	treenum=0;
	treetry=0;
	rootnum=0;
	cycle_count=0;
	nonpolarH=0;
	hflag=0;
	mflag=0;
	Mflag=0;
	aflag=0;
	cflag=0;
	bflag=0;
	rflag=0;
	Aflag=0;
	nanchors = 0;
	nrbonds = 0;
	mindex=0;/*how many lines into mol2 file, bond data starts*/

	xcen=0;
	ycen=0;
	zcen=0;
	newx=0;

	offset=70;
	argindex=1;

	if(argc<4) {
		fprintf(stderr, "Improper syntax for command.  Correct syntax is ...\n");
		fprintf(stderr, "autotors [options] <bond_fname> <input_pdbq_fname> <output_pdbq_fname>\n");
		fprintf(stderr, "or\n");
		fprintf(stderr, "autotors -m [options]  <input_mol2_fname> <output_pdbq_fname>\n");
		fprintf(stderr, "where [options] are...\n");
		fprintf(stderr, "	-m = Use MOL2-formatted file as input\n");
		fprintf(stderr, "	-h = Convert molecule from all-atom to united-atom represenetation\n");
		fprintf(stderr, "	-o = Read old pdbq format file as input (charges in column 55)\n");
		fprintf(stderr, "	-a = Disallow torsions in amide bonds\n");
		fprintf(stderr, "	-b = Disallow torsions in peptide backbone\n");
		fprintf(stderr, "	-c = Add atom connectivity count to ATOM records in output_pdbq_file \n");
		fprintf(stderr, "	-e = <eq_fname> =Use atom type assignments from given file (MOL2 only)\n");
		fprintf(stderr, "	-r = Root automatically set to non-hydrogen atom closest to center of molecule\n");
		fprintf(stderr, "	-M = use  ROTATABLE_BOND and ANCHOR sections to define automatically root and active torsions(valid with MOL2 format only) \n");
		fprintf(stderr, "	-A = Aromatic cycle detection\n");
		exit(1);
	};
	while(argv[argindex][0]=='-'){
		switch(argv[argindex][1]){
			case 'o': offset=55;
				argindex++;
				break;
			case 'h':hflag ++;
				if(hflag>1){
					printf("too many -h flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'm': mflag ++;
				if(mflag>1){
					printf("too many -m flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'M': Mflag ++;
				if(Mflag>1){
					printf("too many -M flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'e': eflag ++;
				if(eflag>1){
					printf("too many -e flags\n");
					exit(1);
				};
				argindex++;
				eqv_file=fopen(argv[argindex++],"r");
				if(eqv_file==NULL){
					printf("Cannot open eqv file\n");
					exit(1);
				}
				tp = eqvtable;
				neqvs = 0;
				while(fgets(line,LINE_LEN,eqv_file)) {
					sscanf(line,"%s%s",tp->sybtyp,tp->autotyp);
					tp++;
					neqvs++;
				}
				fclose(eqv_file);
				break;
			case 'a': aflag ++;
				if(aflag>1){
					printf("too many -a flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'c': cflag ++;
				if(cflag>1){
					printf("too many -c flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'b': bflag ++;
				if(bflag>1){
					printf("too many -b flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'r': rflag ++;
				if(rflag>1){
					printf("too many -r flags\n");
					exit(1);
				};
				argindex++;
				break;
			case 'A': Aflag ++;
				if(Aflag>1){
					printf("too many -A flags\n");
					exit(1);
				};
				argindex++;
				break;
			default:printf("unknown flag in input line\n");
				exit(1);
				break;
		};/*switch on argindex [1]*/
	};
	printf("mol2flag = %d\n",mflag);
	if(mflag) offset=70; 
	printf("hflag = %d\n",hflag);
	if(hflag)printf("offset = %d\n",offset);
	printf("aflag = %d\n",aflag);
	printf("bflag = %d\n",bflag);

/*at this point  need to read in pdbq data into atmptr as well!!!*/
	connect_file=fopen(argv[argindex++],"r");
	if(connect_file==NULL){
		if(mflag){
			printf("Cannot open mol2 file\n");
		} else {
			printf("Cannot open bnd file\n");
		};
		exit(1);
	};
	if(mflag){
		recordtype = UNKNOWNREC;
		while(fgets(line,LINE_LEN,connect_file)) {
			if(line[0] == '@') {
				if(!strncmp(line,"@<TRIPOS>ATOM",13)){	
					recordtype = ATOMREC;
					i = 0;
				}
				else if(!strncmp(line,"@<TRIPOS>BOND",13)){	
					recordtype = BONDREC;
					bp = bndptr;
				}
				else if(!strncmp(line,"@<TRIPOS>ANCHOR",15)){	
					recordtype = ANCHORREC;
				}
				else if(!strncmp(line,"@<TRIPOS>ROTATABLE_BOND",23)){	
					recordtype = RBONDREC;
				}
				else if(!strncmp(line,"@<TRIPOS>MOLECULE",17)){	
					recordtype = MOLECULEREC;
					linenum = 0;
				}
				else recordtype = UNKNOWNREC;
			}
			else {
				switch (recordtype) {
				case MOLECULEREC:
					if(linenum==1) {
						sscanf(line,"%d%d",&natoms,&con_num);
						atmptr=(char *)calloc(natoms+1,LINE_LEN);
						bndptr=(BOND *)calloc(con_num,sizeof(BOND))