diff --git a/onepass.py b/onepass.py
index e7e8617..c0ab217 100644
--- a/onepass.py
+++ b/onepass.py
@@ -1,125 +1,141 @@
 import numpy as np 
 import timeit
 import gurobipy as gp
 from gurobipy import GRB
 
 def addvariables(Z):
     upperbounds=[]
     I=[]
-    for M in range(size_database):
+    for M in database_indices:
         m=len(data['database_ncharges'][M])
         I=I+[(i,j,M) for i in range(m) for j in range(n)]
         upperbounds.append(int(n/m))
 
     x=Z.addVars(I, vtype=GRB.BINARY)
     # dummy variables y associated to number of times molecule M is picked 
-    y=Z.addVars(range(size_database), vtype="I", lb=0, ub=upperbounds)
+    y=Z.addVars(database_indices, vtype="I", lb=0, ub=upperbounds)
     print("Variables added.")
     return x,y,I
 
 def addconstraints(Z,x,y):
     # injection into [n]
     Z.addConstrs(x.sum('*',j,'*') <= 1 for j in range(n))
-    # use at least 60% of target
-    Z.addConstrs(x.sum('*','*','*') >= 0.6*n for j in range(n))
+    # use at least 80% of target
+    Z.addConstr(x.sum() >= 0.8*n)
 
     # additional constraints: take whole molecule or leave it out, and matching charges
-    for M in range(size_database):
+    for M in database_indices:
         CM=data['database_ncharges'][M]
         m=len(CM)
+        # the number of indices of M used (counting multiple use) is at least 50% of its size
+        Z.addConstr(x.sum('*','*',M) >= 0.5*m*y[M])
+        # each index in M is used at most y[M] times
+        Z.addConstrs(x.sum(i,'*',M) <= y[M] for i in range(m))
+        
+        # ignore incompatible charges and atoms of charge 1
         for i in range(m):
-            Z.addConstr(x.sum(i,'*',M) <= y[M])
-            Z.addConstr(x.sum(i,'*',M) >= 0.5*y[M])
             for j in range(n):
                 if(CM[i] != CT[j]):
+                #if(CM[i] != CT[j] or CM[i]==1):
                     Z.addConstr(x[i,j,M]==0)
+        
 
     print("Constraints added.")
     return 0
 
 def setobjective(Z,x,y):
     expr=gp.QuadExpr()
     print("Constructing objective function... ")
-    for M in range(size_database):
+    key=0
+    for M in database_indices:
+        key=key+1
         Mol=data['database_CMs'][M]
         m=len(Mol)
         expr.addTerms(-m, y[M])
         for i in range(m):
             for j in range(m):
                 for k in range(n):
                     for l in range(k,n):
-                        expr.add(x[i,k,M] * x[j,l,M], np.abs(T[k,l]-Mol[i,j])) 
-        print(M+1, "  /  ", size_database)
+                        expr.add(x[i,k,M] * x[j,l,M], np.abs(T[k,l]-Mol[i,j])**2) 
+        print(key, "  /  ", size_database)
 
     Z.setObjective(expr, GRB.MINIMIZE)
     print("Objective function set.")
     return 0
 
 def print_sols(Z, x, y, I):
     SolCount=Z.SolCount
-    for sol in range(SolCount):
+    for solnb in range(SolCount):
         print()
         print("--------------------------------")
-        Z.setParam("SolutionNumber",sol)
-        print("Solution number", sol+1, ", objective value", Z.PoolObjVal)
+        Z.setParam("SolutionNumber",solnb)
+        print("Solution number", solnb+1, ", objective value", Z.PoolObjVal)
         
-        for M in range(size_database):
-            amount_picked=np.rint(y[M].Xn)
-            if y[M].Xn != 0:
-                print("Molecule", data['database_labels'][M], "has been picked", amount_picked, "times ( size", len(data['database_ncharges'][M]),")")
-                m=len(data['database_ncharges'][M])
-                print("Matched indices: ", end='')
+        for M in database_indices:
+            amount_picked=int(np.rint(y[M].Xn))
+            m=len(data['database_ncharges'][M])
+            U=np.zeros((m,n))
+            for k in range(amount_picked):
+                if k==0:
+                    print("Molecule", data['database_labels'][M], "has been picked", amount_picked, "times ( size", len(data['database_ncharges'][M]),")")
+                print(k+1, end=': ')
+                s=0
                 for i in range(m):
                     for j in range(n):
-                        if x[i,j,M].Xn==1:
-                            print(j, end=', ')
-                print()
-
+                        if x[i,j,M].Xn - U[i,j] != 0 and s<0.5*m:
+                            s=s+1
+                            U[i,j]=U[i,j]+1
+                            print(i, "->", j, end=', ')
+                print() 
+        
 def main():
     # construction of the model
     start=timeit.default_timer() 
     Z = gp.Model()
     Z.setParam('OutputFlag',1)
     x,y,I=addvariables(Z)
     addconstraints(Z,x,y)
     setobjective(Z,x,y)
     stop=timeit.default_timer()
     print("Model setup: ", stop-start, "s")
     
     # optimization
     # trying parameters from Z.tune() and others.
     Z.setParam("DegenMoves", 2)
     Z.setParam("ZeroHalfCuts", 2)
     Z.setParam("BQPCuts", 2)
     Z.setParam("RLTCuts", 2)
     
     # time (seconds) and gap (%) limit, (max) number of solutions kept in memory
     Z.setParam("TimeLimit", 180) 
-    Z.setParam("MIPGap", 0.5) # 2-approx. algorithm
+    Z.setParam("MIPGap", 0.5) # 2-approx. algorithm (not really)
     Z.setParam("PoolSolutions", 10)
 
     print("------------")
     print("Optimization")
     print("------------")
     Z.optimize()
 
     print("Optimization runtime: ", Z.RunTime, "s")
     
     if(Z.status == 3):
         print("Model was proven to be infeasible.")
         return 1
     
     print_sols(Z,x,y,I)
     return 0
 
 # global constants 
 
 data=np.load("data.npz", allow_pickle=True)
 target_index=1
 CT=data['target_ncharges'][target_index]
 T=data['target_CMs'][target_index]
 n=len(data['target_ncharges'][target_index])
-size_database=50 # first integer st. model is feasible
-#size_database=len(data['database_ncharges'])
+
+#database_indices=np.loadtxt("filtered"+str(target_index), dtype=int)
+#database_indices=database_indices[:100]
+database_indices=range(100)
+size_database=len(database_indices)
 
 main()