diff --git a/onepass.py b/onepass.py
index 5f976df..e7e8617 100644
--- a/onepass.py
+++ b/onepass.py
@@ -1,121 +1,125 @@
 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):
         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=1)
+    y=Z.addVars(range(size_database), vtype="I", lb=0, ub=upperbounds)
     print("Variables added.")
     return x,y,I
 
 def addconstraints(Z,x,y):
-    # injection into [n], ideally set to equality otherwise the 0 solution is possible
-    Z.addConstrs(x.sum('*',j,'*') == 1 for j in range(n))
+    # 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))
 
     # additional constraints: take whole molecule or leave it out, and matching charges
     for M in range(size_database):
         CM=data['database_ncharges'][M]
         m=len(CM)
         for i in range(m):
-            Z.addConstr(x.sum(i,'*',M)==y[M])
-        for j in range(n):
-            if(CM[i] != CT[j]):
-                Z.addConstr(x[i,j,M]==0)
+            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]):
+                    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):
         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)
 
     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):
         print()
         print("--------------------------------")
         Z.setParam("SolutionNumber",sol)
         print("Solution number", sol+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 i in range(m):
                     for j in range(n):
                         if x[i,j,M].Xn==1:
                             print(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("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=8 # first integer st. model is feasible
+size_database=50 # first integer st. model is feasible
 #size_database=len(data['database_ncharges'])
 
 main()