Proper relaxation method for Si-Ge-Si interface

[DSarkar]

Hi,

I am trying to simulate the interfacial thermal resistance at Si-Ge interfaces, and was following the tutorial on the Si(100) and Si(110) interface using MD. I have attached a file describing the geometry I am currently looking at. From the NPT stabilization of this and similar other previous simulations, it seems to me that the structure is not fully relaxed. In the NPT simulation, there appears to be large destabilizing force at the beginning, which results in weird situations like a wavy character throughout the Si-Ge-Si structure, premature melting at localized regions, etc. In this simulation I set the maximum number of steps in Geometry Optimization as 100000 (maximum possible in VNL), but it ran for 14 fs, and went onto the NPT step. Could you please help me setup the proper relaxation method for this kind of interface? Is there a suitable tutorial or other resource that you can possibly point me to, so that I can get a better hold on the geometry optimization methods? Thanks a lot, I really appreciate your help.

Regards,

Debarghya

[Jess Wellendorff]

I think there are a few things you can try out to improve the calculation:
1) I think you may have misunderstood the purpose of the GeometryOptimization. This is not an MD simulation, it simply minimizes the forces. There is therefore no need to use 100k as maximum number of steps. Try decreasing the maximum force to 0.01 eV/Ang. That relaxation should converge in less than 500 steps.
2) The main issue seems to be the "explosion" you find with MD at 1200 K. I see three different approaches to fix this:
a) Use a pre-equilibration MD run, at e.g.300 K (after the GeometryOptimization).
b) 1200 K may be too much. The tutorial uses 600 K....
c) Your configuration is rather thin along the B direction. This may causes some issues, since phonons couple in all three directions. Try to repeat the configuration along B such that B is roughly equal to A.

[DSarkar]

Hi Jess.

Thanks a lot for your suggestions. One thing that occurred to me was whether the wavy nature of the system and an oscillatory behavior is present because my time step is not sufficiently small to resolve the highest vibrational frequencies of the atoms. In that case, I would want to find the vibrational frequencies and set the time step, say T=0.01*smallest vibrational frequencies. Do you think this is a valid approach? Also, how can this be done in VNL?

Regards,

Debarghya

[Julian Schneider]

I think the time step is not the problem here. You might try 0.5 fs, but even 1 fs might be ok for elements like these.
The wavy nature is probably visually enhanced by the thin dimension of the system, but it doesn't seem unusual at this temperature.

You can get a stable equilibration using the following approach (see attached script):
1. Geometry optimization at constant stress (the off-diagonal lattice components might be slightly non-zero afterwards, they should be reset to zero)
2. Langevin NVT MD at 1200 K
3. NPT simulation at 1200 K

[DSarkar]

Hi Julian,

Thanks a lot for the script and the explanation.

Debarghya