I know it's a boring thing to say, but: It depends on what you want to do.
The way I use Langevin thermostats, is to ensure good equipartitioning in my setup, so that I don't have any local hotspots in the system or something like that. If heat transfer is slow in the system, "global" thermostats like Nosé-Hoover and the like will equilibrate rather slowly. In this context, the damping parameter does not really need to have physical significance, and it only determines how strongly it enforces equipartitioning. That being said, I have found that "1000 time steps" work well for this purpose. For your system, this may translate to a damping constant of 1 ps.
If you wish to model dynamical properties of the system, the issue is much tougher. You correctly note that there exists a connection with the viscosity of the surrounding medium, but this doesn't make all too much sense on/inside a metal. In a liquid, diffusion and viscosity are correlated, but in metals the diffusion of point defects is actually a governed by discrete jumps. So this connection breaks down. Generally speaking I would avoid calculating dynamical properties with Langevin dynamics, and use a global thermostat instead (like
fix nvt or
To return to the first application, there are ways to more systematically optimize the performance of the thermostat, so that it most efficiently samples thermodynamic averages. Whereas a normal Langevin thermostat (
fix langevin or
fix temp/csld in LAMMPS) adds random white noise to the dynamics, you can also add colored noise, using
fix gle. This method implements several ideas developed in the group of Michele Ceriotti. One possible application is "optimal sampling", which guarantees optimal coupling of the thermostat to the system, but there are also options to mimic nuclear quantum effects or vibrational nonequilibrium. There is a dedicated web site that has all the references and documentation.