OK, so you want to measure the impedance of your sensor in a range between 1Hz and 10kHz. From the numbers (10mV, 100p1...100nA) it is between 100kOhm and 100MOhm.
On the data acquisition side, the simplest solution would be to use a 24-bit ADC, so you have enough dynamic range to avoid range switching entirely. If you also want to measure the phase of your impedance accurately, then you'll need a high enough sample rate that the 10kHz sinewave doesn't look like a mess. The DAC that generates the test signal should also have a high enough sample rate to make a clean sine wave. The CPU should also have enough power to process the data.
If you acquire enough samples, then you can also use averaging to reduce noise. Basically, if you acquire a few cycles of sinewave, and multiply this with \$ e^{j \omega t} \$ then average the result, then you can extract phase and magnitude directly. By using averaging, you can extract signal even if it is buried in noise.
If your sensor picks up some ambient 50/60Hz fields, then this method will get rid of them: since it detects the amplitude and phase only at the frequency of interest, it will reject noise and interference at other frequencies. If you simply measure voltage at the transimpedance output, it will not be the case, you will get a lot more noise.
You can also use a log sweep, which is faster.
I got nanoamps resolution with a soundcard, and the current sense resistor was just a few tens of ohms. So the voltage being measured was really tiny, it went through an instrumentation amplifier, then averaged over a bunch of periods. The nice thing about this method is, it gives you a time versus resolution/noise tradeoff. If you find you need lower noise and higher resolution, just use a longer acquisition time.
This would be quite difficult to do with a MSP430, but it is very easy with audio gear. Audio ADC/DAC chips will do the oversampling automatically. In fact, you could use a PC or Raspberry Pi soundcard as data acquisition hardware. Most are AC coupled, so you would need to calibrate out low frequency rolloff or shunt some capacitors, but that's not a problem. Also you probably won't need to code anything, since there are already lots of PC software for impedance analysis... but if you do want to code it yourself, it's easy to do with python.
So the problem is how to turn this tiny current into a voltage suitable for a soundcard input.
Since you're only using AC, the input offset current of your opamp only matters if it's large enough to clip the output. However, your source and feedback impedances are very high, so you should really pay attention to the input noise current of the opamp and its 1/f corner. For example LT1793 has pretty low input current noise.
Note the simulation trace in the question shows the opamp is clipping. If you use a single supply, then its input common mode and output voltage range need to include ground.
It would probably be easier to use positive and negative supplies, to avoid clipping the opamp at zero input current if you get a specimen with an offset voltage that goes in the wrong direction. Also if you use AC current to measure your impedance, then it is much easier to use split supplies for the opamp. Otherwise you'll have to AC-couple with a capacitor to a reference voltage somewhere around midsupply, and then think about the cap's leakage current, etc.
The other answers provide enough info about how to do the TIA part, so I won't elaborate.
Edit: DSP stuff
To plot modulus and phase of impedance vs frequency, you can use two methods:
This is the classic test, it sweeps the entire frequency range of interest. Detailed explanations and math are in the link, it's about 1-2 pages of python code. Since it requires FFT of the entire waveform, it has to fit in memory, so this isn't possible with a microcontroller. It is the fastest method if you are interested in the whole range of frequencies.
- Run multiple single-frequency tests, stepping frequency
This is simpler. If you're only interested in a few frequency points, then this is a good option. It can have lower noise than the swept sine if it spends more time averaging one frequency, but of course if you're interested in the full frequency range, it will be much slower since each frequency has to be measured independently, and that includes a bit of settling time after stepping to the next frequency.
It works like this:
Generate sinewave \$ Ao sin(\omega t) \$ and play it with DAC
Record measurement on device under test, it'll also be a sinewave with different amplitude and phase shift. Say the recorded data is \$ Ar sin(\omega t + \phi) \$ with Ar=recorded amplitude and \$ \phi \$ phase shift relative to the original sinewave.
So you compute \$ e^{j \omega t} = cos \omega t + j sin \omega t\$, multiply that with the received signal, and average the result over an integer number of sine periods. The sine period doesn't have to be an integer number of samples, but the averaging should be over an integer number of periods.
Multiplying these two and averaging is just like calculating the discrete time integral of the product of recorded signal with \$ e^{j \omega t} \$. If you write down the integral and solve it, you'll notice the result is a complex number whose amplitude is Ar and phase is \$ \phi \$.
Or maybe it's \$ -\phi \$, I don't remember, but you get the idea.
That's basically how a RF detector works, it multiplies the incoming wave with sin and cos, then averages.
So you get a complex number \$ A(f) \$ that encodes amplitude and phase for each frequency. The swept sine also gives results in the form of complex numbers, so from now on, the calibration method is the same. You have to keep the result as complex number, because that will be useful.
In order to calibrate it, you have to focus on the fact that what you want to measure is the difference between the electrochemical sensor working in the chemical solution, and an identical sensor that is not dipped in the solution. Or maybe it's something else, that's for you to decide, but the point of calibration is to remove all the stuff you don't want to be part of your measurement.
For example you have to cancel all the phase shift due to the setup, DACs, ADCs, sampling delays, buffering, etc (which will not depend on frequency) and the phase shift due to capacitance and analog filters (which will depend on frequency). And you also have to measure a known resistor to know what an ADC amplitude means in terms of analog values.
Since you're measuring a resistor up to 100 Mohms at up to 10kHz, and a parasitic capacitance of 1pF will have 15Mohms impedance at 10kHz, you'll have to be very careful about capacitance, otherwise it will completely swamp your measurement. It's not possible to get rid of capacitance, but if you make it constant, then you can calibrate it out. If picofarads matter, this means you can't have flying wires that could be set in different shapes or with different spacings, because that would change the capacitance. You'll have to use coax. Or you could use no cables at all, and put the TIA on a small PCB right next to the sensor, for example.
Another reason to not use multiple ranges on your ADC is that your setup will probably have like at least 10pF in parallel with the sensor, so at 10kHz on the 100pA range the current through the capacitance will be a lot higher than what you try to measure, so if you use range switching it will just clip the opamp and ADC. If you use an ADC with lots of bits, then no problem.
If you manage to make the sensor capacitance constant and repeatable, then you can fiddle with the TIA feedback cap to somewhat compensate for it. It won't be perfect, unless you manage to build all the sensors absolutely identical, or you make the TIA part of the sensor assembly so each sensor can have its personalized feedback cap, but it should help.
To calibrate you need to be able to replace the sensor with a resistor, so you'll need a connector.
So first you remove the sensor and measure an open circuit, this will give you the response of your setup including all parasitic capacitance, so for this to be accurate, all the wiring must be the same, just the sensor removed. Or you could use a dummy electrochemical sensor without liquid so it has infinite resistance but it still has the capacitance it would have when used in the chemical solution.
Then you can replace the sensor with a known resistor value and measure again. Or, for better results, you would wire a known resistor value in parallel with a sensor that is not dipped in solution, so the only change is resistance, not capacitance.
The method is similar to 2-port network analyzer calibration, but you'll probably be able to ignore most of the terms and use only two calibration points, open and "known load". You can look for literature on that, or you can draw schematics of the various combinations of sensors, no sensors, and calibration resistors... then calculate the impedance of that, put an "=" sign and the measured data on the other side, and solve that system of equations.
If you get the math right, the end result is just a bunch of complex multiplications, divisions and substractions between measured data and calibration data, so it's not difficult to program.