anchor.diff:  contains
     
     subroutine rattle_pos
  -  double precision :: g !correction factor 
  +  !double precision :: g !correction factor 
     double precision :: s(3) !direction vector 
     double precision :: r(3) !old direction vector
     double precision :: d
  -  
  +  !double precision :: q_1(3),q_2(3)  
  +
  +
     d  = 2.d0*sigma + 0.5d0 !distance between the colloids in the dimer
   
     s = colloids% pos_old(:,1) - colloids% pos_old(:,2)
     r = pos_rattle(:,1) - pos_rattle(:,2)
     g = (dot_product(s,s) - d**2)/(2.d0*(dot_product(s,r)*(2.d0/mass)))
     
  +  !q_1 = colloids% vel(:,1) + colloids% force(:,1)*dt/2/mass
  +  !q_2 = colloids% vel(:,2) + colloids% force(:,2)*dt/2/mass
  +
  +  !g = dt*(-(2*dt**2*dot_product(q_2,r)/mass - 2*dt**2*dot_product(q_1,r)/mass + 2*dt*dot_product(q_2,r)/mass &
  +  !    - 2*dt*dot_product(q_1,r)/mass)) &
  +  !   - sqrt((2*dt**2*dot_product(q_2,r)/mass - 2*dt**2*dot_product(q_1,r)/mass + 2*dt*dot_product(q_2,r)/mass &
  +  !    - 2*dt*dot_product(q_1,r)/mass)**2 &
  +  !    -4*dt**2*(dot_product(r,r)*2/mass**2)*(dot_product(colloids% pos(:,1),colloids% pos(:,1)) &
  +  !    + dot_product(colloids% pos(:,2),colloids% pos(:,2)) - d**2 &
  +  !    + dt**2*(dot_product(q_1,q_1) + dot_product(q_2,q_2)) &
  +  !    + 2*dt*(dot_product(colloids% pos(:,1),q_1)+dot_product(colloids% pos(:,2),q_2))) &
  +  !    /(dt**2*(dot_product(r,r)*2/mass**2)))
  +
     colloids% pos_old(:,1) = colloids% pos_old(:,1) - g*r/mass
     colloids% pos_old(:,2) = colloids% pos_old(:,2) + g*r/mass
     
  
positive.code:   integer :: j, k
    integer, allocatable :: seed(:)
  
    call random_seed(size = seed_size)
    allocate(seed(seed_size))
    call system_clock(count=clock)
  
negative.code: 
          call switch(solvent% force, solvent% force_old)
          call switch(colloids% force, colloids% force_old)
  
          solvent% force = 0
          colloids% force = 0
          e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
  