Hi Jenny,
Unless you are explicitly interested in thermal or temperature dependent properties (e.g. thermal expansion), performing a geometry optimization at 0K and then calculating the desired properties based on the final structure is in fact the most common approach, in particular when you are interested in electronic-structure-related properties.
I see that your MD temperature is 50 K, and for this low temperature I would expect the difference to the 0K-structure to be very small.
The disadvantage of finite-temperature-structures is that there is no unique structure but rather an ensemble of configurations that you would have to average over.
If you really want to use finite temperature structures, then in your case, you would have to use the the NPTMelchionna-thermostat isntead of the NPTBerendsen-thermostat. The reason is that you have an anisotropic structure that is only periodic in z-direction, and NPTBerendsen treats everything as isotropic (cf. MD-tutorial), which probably causes your error message.
So you'd have to use something like:
method = NPTMelchionna(
time_step=1*femtoSecond,
reservoir_temperature=50*Kelvin,
external_stress=1*bar,
thermostat_timescale=100*femtoSecond,
barostat_timescale=100*femtoSecond,
bulk_modulus=1e+06*bar,
initial_velocity=initial_velocity,
mask=[[False, False, False], [False, False, False], [False, False, True]]
)
