anchor.diff:    dt = tau / N_MD_steps
     N_loop = PTread_i(config, 'N_loop')
     N_therm = PTread_i(config, 'N_therm')
  +  vxz_interval = PTread(config, 'vxz_interval')
   
     wall_t = PTread_dvec(config, 'wall_T', 2)
     T = PTread_d(config, 'T')
  
positive.code: 
    end subroutine rescale_at_walls
  
    subroutine compute_vxz
      integer :: i, s, ix, iy, iz
  
      do i = 1, solvent%Nmax
  
negative.code:         rhoz% data = rhoz% data / rhoz% dx
          call elem_rhoz% append(rhoz% data, i, i*tau)
          rhoz% data = 0
          call compute_vxz_and_vxyz
          if (modulo(i-N_therm, vxz_interval)==0) then
             call div_vxz()
             call v_xz_el%append(v_xz, i, i*tau)
             v_xz_count = 0
             v_xz = 0
  