Calculating force constants for very large systems

[michael_walkup]

I'm trying to calculate the force constants for a system of 324 ions so that I can compute the phonon dispersion and DOS. I've done some test runs using 4 nodes and 32 tasks per node with a 2x2x2 k-point and q-point mesh, and in four hours the calculation only makes it through 79/972 degrees of freedom. I could increase the computing time and nodes requested, but I'm wondering if there are some optimizations I can utilize to reduce the resources needed. What is the best way to optimize this calculation? Would using DFPT instead of finite differences be more efficient? I have attached the INCAR, POSCAR, POTCAR, KPOINTS, QPOINTS, output and batch script files for reference.

-Michael

[christopher_sheldon1]

Hi Michael, I'm afraid that I don't see your files attached. Could you reattach them?

I doubt that DFPT would be more efficient. Have you already relaxed the structure or is this what you are performing?

[michael_walkup]

My apologies, I could've sworn I attached the zip folder in the last message. The structure is already relaxed, yes.

Edit: Upon trying to attach it again, the attachment status indicates:
"ERROR
Sorry, the board attachment quota has been reached."
So I'm not sure what to do in this instance.

-Michael

[christopher_sheldon1]

It looks like an error our end. We'll work out a solution and then get back to you.

Chris

[christopher_sheldon1]

We've removed some old files so there should be more space now. Could you try re-uploading?

Chris

[michael_walkup]

Looks like it's working now. Thanks!

-Michael

[christopher_sheldon1]

Hi Michael,

Thanks for uploading them! Glad we've got it working

Your system is a sort of dibutoxyindole-based porphyrin crystal. I'm presuming that this is the relaxed (or experimental) structure? If you could upload your OUTCAR file, this might add clarity for a few points. The first is, since it is an organic crystal, I presume that it is an insulator, rather than a conductor?

A few minor points to consider, which are essentially optional:
If it is an insulator, then you can change the smearing from (ISMEAR = 2, SIMGA = 0.2) to (ISMEAR = 0, SIGMA = 0.05), i.e. 2nd order Methfessel-Paxton to Gaussian. This is a fairly minor point though.
A second, is that you are using IBRION = 6 but have switched symmetry off except for sampling of the Brillouin zone (ISYM = 1). You could try turning symmetry on as your system looks pretty symmetric.
It might also be worth adding in a disperison correction (e.g. D2, D3, cf. IVDW), as for an organic crystal this tends to be quite important.

Then a few bigger ones that might make quite a large difference to your results:
Your POTIM is quite large (0.8 ), I would change this to be 0.015, although the software may do this automatically already.
A big one is that you have set your convergence criteria for the forces to be extremely tight (EDIFFG = -1E-5 eV/A) - the minus indicates forces, rather than energy difference is the criterion. I would loosen this significantly to (EDIFFG = -0.01).

In terms of time-saving, there are a few suggestions that I have:
Firstly, you could reduce the k-point mesh from 2x2x2 to just the gamma point. This would offer a significant time-saving, although you would have to test how it changes the phonons/DOS.
Secondly, you could change PREC from Accurate to Normal. This should give a big time-saving.
Finally, I noticed that your structure contains 2 porphyrin rings. It might be possible to find a symmetrically equivalent structure with only one in. That would save a lot of time. However, this depends a lot on your system.

Does this answer your question?

Best wishes,

Chris

[michael_walkup]

Thank you very much for your help, Chris. Yes, the crystal is an insulator, so I will include your tweaks to ISMEAR. I've attached the OUTCAR file if that helps clear things up. The advice about using just the gamma point seems especially useful; I had just assumed that a 2x2x2 mesh was the minimum. I'll look into seeing if modifying the structure itself would be feasible, but first I'm going to try your suggestions and see if that may be sufficient.

-Michael

[christopher_sheldon1]

You're welcome. Let me know how the calculation goes.

Thanks for the OUTCAR file (I've downloaded but have removed it here to ensure that there's enough space on the forum to upload). It definitely looks like an insulator, although it does have an interesting structure with a doubly degenerate occupied band, which I presume is from weakly bound/ non-bonding pi-electrons in the aromatic ring. Although that is just speculation.

Best wishes,

Chris

[christopher_sheldon1]

Hi Michael,

I just had a conversation about symmetry with a colleague and your question came to mind. You can test the symmetry by looking at the symmetry and precision section for the POSCAR tag on the wiki. There's a link to a script there that can test this and give you a primitive unit cell, if available.

Additionally, if you do the gamma point, rather than a 2x2x2 k-point mesh, you can save additional time by using the gamma point-only version of VASP, vasp_gam, rather than the standard vasp_std. This should give you a further time saving.

Best,

Chris