How about if you just add 5e-9 to all numbers:
close all
gm = 1e-6:1e-3:1;
grid on
hold on
%Rd = 1K
einv = 2.*sqrt((4e-9)^2+(4.*1.38e-23.*298.15.*1e3))
einv = einv./(gm*1e3) + 5e-9;
plot(gm,einv,'color',[0,0,0], 'LineWidth', 2);
%Rd = 10K
einv = 2.*sqrt((4e-9)^2+(4.*1.38e-23*298.15*10e3))
einv = einv./(gm*10e3) + 5e-9;
plot(gm,einv,'color',[0.25,0.25,0.25], 'LineWidth', 2);
%Rd = 100K
einv = 2.*sqrt((4e-9)^2+(4.*1.38e-23*298.15*100e3))
einv = einv./(gm*100e3) + 5e-9;
plot(gm,einv,'color',[0.5,0.5,0.5], 'LineWidth', 2);
%Rd = 1M
einv = 2.*sqrt((4e-9)^2+(4.*1.38e-23*298.15*1e6))
einv = einv./(gm*1e6) + 5e-9;
plot(gm,einv,'color',[0.75,0.75,0.75], 'LineWidth', 2);
xlabel ('gm (S)'); ylabel ( 'EINV_{min} (V/sqrt(Hz))');
set(gca,'xscale','log')
set(gca,'yscale','log')
legend('1K\Omega','10K\Omega','100K\Omega','1M\Omega');
%print('einv','-dpng')
Looks fairly good to me.
