diff --git a/bwape.c b/bwape.c
index 63783c8..204bbfe 100644
--- a/bwape.c
+++ b/bwape.c
@@ -395,16 +395,19 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
 		if (opt->N_multi || opt->n_multi) {
 			for (j = 0; j < 2; ++j) {
 				if (p[j]->type != BWA_TYPE_NO_MATCH) {
-					int k;
+					int k, n_multi;
 					if (!(p[j]->extra_flag&SAM_FPP) && p[1-j]->type != BWA_TYPE_NO_MATCH) {
 						bwa_aln2seq_core(d->aln[j].n, d->aln[j].a, p[j], 0, p[j]->c1+p[j]->c2-1 > opt->N_multi? opt->n_multi : opt->N_multi);
 					} else bwa_aln2seq_core(d->aln[j].n, d->aln[j].a, p[j], 0, opt->n_multi);
-					for (k = 0; k < p[j]->n_multi; ++k) {
+					for (k = 0, n_multi = 0; k < p[j]->n_multi; ++k) {
 						int strand;
 						bwt_multi1_t *q = p[j]->multi + k;
 						q->pos = bwa_sa2pos(bns, bwt, q->pos, p[j]->len, &strand);
 						q->strand = strand;
+						if (q->pos != p[j]->pos)
+							p[j]->multi[n_multi++] = *q;
 					}
+					p[j]->n_multi = n_multi;
 				}
 			}
 		}
diff --git a/bwase.c b/bwase.c
index 76480d4..8fa8d45 100644
--- a/bwase.c
+++ b/bwase.c
@@ -84,12 +84,6 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
 			}
 		}
 		s->n_multi = z;
-		/*// the following code removes the primary hit, but this leads to a bug in the PE mode
-		for (k = z = 0; k < s->n_multi; ++k)
-			if (s->multi[k].pos != s->sa)
-				s->multi[z++] = s->multi[k];
-		s->n_multi = z < n_multi? z : n_multi;
-		*/
 	}
 }
 
@@ -141,19 +135,23 @@ void bwa_cal_pac_pos_core(const bntseq_t *bns, const bwt_t *bwt, bwa_seq_t *seq,
 
 void bwa_cal_pac_pos(const bntseq_t *bns, const char *prefix, int n_seqs, bwa_seq_t *seqs, int max_mm, float fnr)
 {
-	int i, j, strand;
+	int i, j, strand, n_multi;
 	char str[1024];
 	bwt_t *bwt;
 	// load forward SA
 	strcpy(str, prefix); strcat(str, ".bwt");  bwt = bwt_restore_bwt(str);
 	strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt);
 	for (i = 0; i != n_seqs; ++i) {
-		bwa_cal_pac_pos_core(bns, bwt, &seqs[i], max_mm, fnr);
-		for (j = 0; j < seqs[i].n_multi; ++j) {
-			bwt_multi1_t *p = seqs[i].multi + j;
-			p->pos = bwa_sa2pos(bns, bwt, p->pos, seqs[i].len, &strand);
-			p->strand = strand;
+		bwa_seq_t *p = &seqs[i];
+		bwa_cal_pac_pos_core(bns, bwt, p, max_mm, fnr);
+		for (j = n_multi = 0; j < p->n_multi; ++j) {
+			bwt_multi1_t *q = p->multi + j;
+			q->pos = bwa_sa2pos(bns, bwt, q->pos, p->len, &strand);
+			q->strand = strand;
+			if (q->pos != p->pos)
+				p->multi[n_multi++] = *q;
 		}
+		p->n_multi = n_multi;
 	}
 	bwt_destroy(bwt);
 }
diff --git a/main.c b/main.c
index d15d14f..4d4402f 100644
--- a/main.c
+++ b/main.c
@@ -4,7 +4,7 @@
 #include "utils.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.6.0-r89-dev"
+#define PACKAGE_VERSION "0.6.0-r90-dev"
 #endif
 
 static int usage()