2.1 PhonoPy to PyEPFD conversion
Let us asuume the the root path of the pyepfd is pyepfd_path/ . You will find pyepfd_path/utils/phonopy/phonopy2pyepfd.py file to convert FORCE_CONSTANTS file into PyEPFD format. Let us first try to run it.
[1]:
%%bash
#Below I'm defining my pyepfd_path relative to tutorials directory
pyepfd_path=../../../
#Now I am moving to the directory which contains FORCE_CONSTANTS file;
#The phonopy2pyepfd script must be run from the directory which contains FORCE_CONSTANTS file
cd PhonoPy_DynMatrix/
python3 $pyepfd_path/utils/phonopy/phonopy2pyepfd.py
cd ../
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2024-12-21 14:28:05.715950
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
Usage: python phonopy2pyepfd.py path_to_opt_geometry.xyz path_to_pyepfd_restart.xml <freq_scale_factor>[optional]
We have not provided any command line arguments. Therefore the code is complaining and asking for command line arguments. (1) The first argument is the full path of the xyz file containing optimized geometry. (2) The second argument is the full path of the pyepfd_restart file you want to create. (3) Third argument is optional and is necessary if you want to scale the frequency as explained below. Here we computed the dynamical matrices and frequencies of the NV-center in diamond using PBE functional. Then for a much smaller supercell, the frequencies of pristine diamond is calculated using both DDH and PBE. The ratio of highest frequency as obtained from DDH and PBE is found to be 1.04. Now if all PBE frequencies are just scaled by 1.04, it is found the VDOS of diamond is well reproduced (see Yu Jin et al, Phys. Rev. Mater. 2021, 5, 084603). We used this approach to estimate the DDH frequencies and therefore as a third argument we will use 1.04. Lets convert Phonopy FORCE_CONSTANTS to pyepfd_restart file which we name: nv216_ddh_phonon.xml
[2]:
%%bash
#Below I'm defining my pyepfd_path
pyepfd_path=/home/arpan/Work/Code-developments/EPFD/pyepfd-dev
#Now I am moving to the directory which contains FORCE_CONSTANTS file;
#The phonopy2pyepfd script must be run from the directory which contains FORCE_CONSTANTS file
cd PhonoPy_DynMatrix/
python3 $pyepfd_path/utils/phonopy/phonopy2pyepfd.py NV216-DDH-opt.xyz nv216_ddh_phonon.xml 1.04
cd ../
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2024-12-21 14:28:26.756239
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
Time spent on xyz class: 0.0005693435668945312 s.
#Mode Frequency(cm-1)
1 -2.012360757094569e-05
2 -1.4041660928339495e-05
3 -1.3514296592496058e-05
4 338.9944731789015
5 338.99447317890167
6 341.32347390167575
7 342.0295327556809
8 342.92474739291674
9 342.92474739291686
10 346.400452538121
11 346.40045253812167
12 346.48592041988337
13 350.0973196607621
14 351.5584805116746
15 351.558480511675
16 420.87626190615435
17 420.8762619061546
18 440.81765239916786
19 442.423066448968
20 442.4230664489681
21 442.67618191164087
22 444.68584051871096
23 445.0341807809395
24 445.61838187185685
25 445.61838187185685
26 447.1099646624555
27 447.1099646624561
28 465.76194706301123
29 465.7619470630116
30 466.7288106338762
31 468.12460633396955
32 468.1246063339696
33 470.87322719459223
34 470.87322719459223
35 471.9826611263647
36 472.1764574680828
37 473.74096511180477
38 473.7409651118056
39 474.06574111224785
40 477.07811385876516
41 477.07811385876533
42 478.20056310581975
43 478.20056310582015
44 499.4655965516689
45 503.5106787725517
46 503.5106787725517
47 507.7794995565502
48 508.96299922590305
49 508.9629992259032
50 513.1426019300271
51 513.640132795036
52 513.6401327950366
53 525.8211422511687
54 525.8211422511687
55 526.1293210973066
56 527.5956468104731
57 527.5956468104735
58 527.9280967809112
59 529.9113742835578
60 529.9113742835582
61 534.6260190061872
62 586.17348704413
63 592.4593681927196
64 592.4593681927198
65 605.1873190273433
66 605.1873190273437
67 608.5375881747462
68 612.6449250126484
69 613.6542452610406
70 613.654245261041
71 617.2239747403517
72 617.353012711691
73 618.1433265109004
74 618.1433265109005
75 618.9625200304191
76 619.3436179779969
77 619.3436179779974
78 619.8531525612515
79 619.8531525612517
80 620.3778003699344
81 621.6720787931417
82 622.273431364476
83 622.2734313644763
84 622.5914547120494
85 622.637877613136
86 622.6378776131362
87 641.1619133733384
88 641.1619133733386
89 642.4346548319653
90 654.0554127129856
91 654.0554127129858
92 654.7550444327242
93 657.5853895971659
94 658.3636360910316
95 658.3636360910317
96 659.1154496493292
97 659.1154496493294
98 659.3142579601241
99 667.3598637404958
100 667.3598637404959
101 669.7509015373637
102 669.7509015373638
103 670.7650705314965
104 675.7278190482306
105 675.7278190482308
106 675.8387132152956
107 679.8518703244881
108 679.8518703244881
109 683.1065036664543
110 683.4912864303413
111 683.4912864303415
112 685.8433160456997
113 685.8433160456999
114 687.8645500341812
115 695.0383667840517
116 697.8495881226048
117 697.849588122605
118 698.3475021184203
119 700.8618968715775
120 700.8618968715776
121 701.9584347263983
122 702.0975268710134
123 702.473871492897
124 702.473871492897
125 702.6379265425905
126 702.6379265425907
127 702.8428066098527
128 703.6767274774601
129 703.6767274774604
130 703.8480109027374
131 703.8854327463605
132 703.9618183116603
133 703.9618183116604
134 704.7675971838077
135 705.1286647994029
136 705.1286647994029
137 706.454229997376
138 706.4542299973764
139 706.784371868036
140 708.6907986852392
141 708.6907986852393
142 709.2355334279293
143 709.2355334279301
144 709.3864034742547
145 710.2893365388633
146 711.3215401896007
147 711.5119309425074
148 711.7118087725946
149 711.7118087725948
150 712.6789867279547
151 712.6789867279548
152 715.5920032900502
153 715.6610014235569
154 715.6610014235569
155 717.6710446427047
156 717.6710446427052
157 720.3768680248771
158 720.9074008697233
159 720.9903391248296
160 720.99033912483
161 722.9587402586866
162 722.958740258687
163 723.9645149837437
164 723.9668861105919
165 723.966886110592
166 726.3153605219288
167 727.6954888789022
168 727.6954888789023
169 729.9559846038039
170 757.5633342557707
171 761.5292112478344
172 761.5292112478345
173 763.9940899236308
174 773.726908605764
175 773.7618996128731
176 773.7618996128731
177 775.2966826675402
178 775.3144767388245
179 775.3144767388246
180 775.5490828225627
181 775.549082822563
182 777.2144153869168
183 777.5273440572663
184 777.5273440572663
185 778.4370814385526
186 779.8831785942879
187 779.883178594288
188 780.5695605738811
189 781.714001358905
190 782.2933390176347
191 782.2933390176349
192 782.8621428758133
193 782.8621428758134
194 792.3187778006892
195 792.6174434689083
196 792.6174434689084
197 793.5657195930971
198 793.5657195930976
199 795.100084161615
200 809.699879769142
201 858.6387970627114
202 858.6387970627117
203 858.6709676309216
204 873.2408737943707
205 873.2408737943709
206 875.5194571324957
207 875.5194571324961
208 876.2582262177588
209 877.786135635478
210 878.8315370134864
211 878.8315370134864
212 879.1939897802633
213 880.0379241288211
214 880.0453411242968
215 880.0453411242968
216 880.4716510811276
217 880.9164432473942
218 880.9172150762377
219 880.9172150762379
220 882.994148117532
221 882.9941481175322
222 884.0648627103276
223 884.2165087906973
224 884.6264537457316
225 884.626453745732
226 885.520632051517
227 885.5206320515174
228 885.6933099730941
229 886.340308216929
230 886.8844635758802
231 886.8844635758802
232 891.0185441786888
233 909.9278623254145
234 910.7293970402704
235 910.7293970402704
236 911.970259390193
237 911.9702593901931
238 912.7304047053342
239 913.2615989106265
240 913.5580140432611
241 913.5721256366639
242 913.5721256366642
243 914.9034687396849
244 914.9034687396849
245 915.2223007822811
246 915.325767550889
247 915.8225829118479
248 915.8225829118479
249 916.2541618710834
250 916.2541618710834
251 917.1724797028711
252 917.1724797028713
253 920.3579267866447
254 924.2836950277767
255 924.2836950277768
256 945.4012549482727
257 951.2804955987233
258 951.2804955987233
259 956.6395423659086
260 958.3518241542702
261 958.3518241542703
262 963.8573481378633
263 976.9882279305699
264 976.9882279305699
265 978.5586042386153
266 980.3173466533042
267 980.3173466533043
268 982.5650295652307
269 982.5650295652309
270 982.719436050856
271 982.7520181844908
272 983.3686561225842
273 983.3686561225842
274 983.6520592368813
275 983.6520592368813
276 984.236028783244
277 984.236028783244
278 984.9799843171774
279 986.6973742206962
280 987.475343747929
281 989.3220786773351
282 989.7083753993882
283 989.7083753993883
284 990.3679415689517
285 990.3679415689519
286 1006.4645227982145
287 1013.7332026487442
288 1013.7332026487442
289 1014.9331796050901
290 1016.2181321069562
291 1016.6533238790976
292 1016.6533238790978
293 1017.4145387016151
294 1017.4145387016151
295 1019.4726259995647
296 1019.4726259995648
297 1019.6550014129574
298 1020.3748407011531
299 1020.3748407011531
300 1021.0576933626946
301 1021.0576933626946
302 1022.501506787343
303 1022.5246180679937
304 1022.5246180679943
305 1022.9689502483919
306 1024.4396174722476
307 1025.8134349258416
308 1026.873965682353
309 1026.873965682353
310 1036.827657284932
311 1038.7306513629803
312 1038.7306513629803
313 1042.0166668649972
314 1042.0166668649972
315 1042.4854169845157
316 1045.2843422656813
317 1046.2630201768688
318 1046.4759220445942
319 1046.4759220445942
320 1047.668749717824
321 1047.668749717824
322 1048.1514433957393
323 1050.4215138010454
324 1051.3819500456748
325 1051.3819500456748
326 1052.959586772436
327 1054.0118925447061
328 1054.0118925447064
329 1055.0600115978393
330 1055.0600115978395
331 1055.6858980332868
332 1055.685898033287
333 1057.5765127720074
334 1078.1648400910228
335 1078.164840091023
336 1078.88619622168
337 1081.955637159362
338 1084.0198136316367
339 1084.0198136316376
340 1085.864056352729
341 1091.0592545255163
342 1091.0592545255165
343 1091.8832011092159
344 1092.0888272019936
345 1092.0888272019936
346 1094.2938062063556
347 1094.2938062063558
348 1094.9045646924858
349 1095.265221902135
350 1096.889828996235
351 1096.8898289962353
352 1116.28374728688
353 1120.1055510992587
354 1120.105551099259
355 1121.4383947842937
356 1121.4383947842937
357 1128.3174754832457
358 1130.5300174909885
359 1130.5300174909892
360 1131.5457655773757
361 1131.9076021499027
362 1131.907602149903
363 1132.2013869956766
364 1141.040505849564
365 1141.0405058495642
366 1144.1189038747511
367 1144.8459700656526
368 1144.8459700656529
369 1146.802983083695
370 1148.683982397188
371 1148.6839823971882
372 1151.5464823733687
373 1153.1472005017642
374 1153.6230491980095
375 1153.6230491980098
376 1154.381802905956
377 1155.2319447993573
378 1155.2319447993573
379 1155.4245951033317
380 1156.191401783216
381 1156.191401783216
382 1158.4615761382943
383 1158.5603214714579
384 1158.5603214714579
385 1159.7863634406572
386 1159.7863634406574
387 1159.866592582785
388 1160.5619541620952
389 1161.8788485477091
390 1161.8788485477091
391 1162.7175142764188
392 1167.545962638237
393 1167.545962638237
394 1168.2263869403264
395 1171.3512164975623
396 1171.3512164975623
397 1172.6878760019622
398 1173.0333226719104
399 1173.0333226719106
400 1175.011457057081
401 1178.4779938092859
402 1178.477993809286
403 1181.529122339759
404 1183.5806279094113
405 1183.8654795075317
406 1183.8654795075317
407 1184.519759631475
408 1184.519759631475
409 1187.8688633787235
410 1189.2087796716394
411 1189.2087796716394
412 1189.7081362216131
413 1190.3275734078277
414 1190.327573407828
415 1190.4373458969417
416 1190.6721285840308
417 1190.6721285840308
418 1191.0201073719822
419 1191.2314070061964
420 1191.784865183616
421 1191.784865183616
422 1192.4830470214515
423 1192.7223162319997
424 1192.722316232
425 1193.2307784884263
426 1193.2307784884263
427 1193.5258873349758
428 1194.112783815516
429 1194.343612575905
430 1194.343612575905
431 1194.61425494013
432 1194.9623479228296
433 1194.9623479228298
434 1195.4691352630352
435 1195.4691352630355
436 1195.544519763195
437 1195.6460220868273
438 1195.8939458318227
439 1195.8939458318227
440 1197.0152661055522
441 1197.0152661055522
442 1197.7948219853222
443 1202.8988549714645
444 1202.9573544325906
445 1202.9573544325906
446 1204.854679135154
447 1205.3143407453003
448 1205.3143407453006
449 1206.3065194364126
450 1206.3065194364126
451 1208.2962815208782
452 1209.6816351420894
453 1209.6816351420894
454 1211.5472598362207
455 1211.547259836221
456 1213.9372057110734
457 1214.2859115308222
458 1214.8183796532242
459 1214.8183796532242
460 1216.2095597073403
461 1216.2095597073403
462 1216.704597940631
463 1216.9181720253816
464 1216.9181720253816
465 1217.5641740466056
466 1217.8649913839795
467 1217.8649913839797
468 1218.8442197611053
469 1218.8442197611053
470 1218.9145042942491
471 1219.5073942244935
472 1220.8337793289834
473 1220.8337793289834
474 1221.019368853806
475 1222.5890149209906
476 1222.5890149209906
477 1222.8303117936186
478 1223.5402755716925
479 1224.4416760003023
480 1224.778377336408
481 1224.7783773364085
482 1225.2325038548595
483 1225.2325038548595
484 1226.0681963348711
485 1226.0681963348713
486 1226.7777110984093
487 1226.9705021851914
488 1228.0618640227021
489 1228.30035142026
490 1228.3003514202603
491 1228.7219177668014
492 1228.7219177668014
493 1230.6999439176518
494 1231.164296709747
495 1231.44602181188
496 1231.4460218118804
497 1234.1143746723333
498 1234.1143746723333
499 1236.6875899313602
500 1236.6875899313604
501 1237.1467488909252
502 1244.841459866062
503 1246.9383426667396
504 1246.9383426667396
505 1247.8947923659991
506 1247.8947923659991
507 1248.4568547247138
508 1248.8100410856262
509 1248.8100410856262
510 1249.074208078197
511 1251.517599870529
512 1251.7524184154634
513 1251.7524184154634
514 1252.813296062043
515 1252.813296062043
516 1252.9662153973995
517 1253.2207971355708
518 1253.2207971355713
519 1254.0005919500286
520 1254.2495819980284
521 1254.4958759859683
522 1254.4958759859683
523 1255.1691091410694
524 1255.1691091410696
525 1256.1248294128538
526 1256.3160309084735
527 1256.3160309084735
528 1257.4878195792232
529 1258.0333261407677
530 1258.0333261407677
531 1259.0025064644105
532 1259.7484848790473
533 1259.7484848790473
534 1261.2010286226935
535 1261.201028622694
536 1261.9498256597087
537 1262.2198292120238
538 1262.2198292120243
539 1262.3392739367991
540 1262.3392739367998
541 1262.5987410607977
542 1262.607725205875
543 1262.7168148951669
544 1262.716814895167
545 1263.3287130984806
546 1263.3287130984809
547 1263.8226629733474
548 1263.8226629733476
549 1264.17675040753
550 1264.17675040753
551 1264.664967097333
552 1264.9483540714868
553 1265.6228133275827
554 1265.9953629868392
555 1265.9953629868392
556 1266.4493686510216
557 1266.4493686510218
558 1266.5888049550406
559 1266.7527157096486
560 1267.0680407734674
561 1267.436652598392
562 1267.4366525983921
563 1267.6231531515978
564 1267.623153151598
565 1268.3391256931186
566 1268.453264692743
567 1268.4540861538314
568 1268.4540861538314
569 1268.775038468763
570 1269.113750131574
571 1269.579788395582
572 1269.5797883955822
573 1269.9979914085973
574 1270.5156086129518
575 1270.515608612952
576 1270.8743569689666
577 1271.8196322979986
578 1271.8196322979989
579 1273.1829664078848
580 1273.1829664078853
581 1275.287762622339
582 1276.5038804827982
583 1276.5038804827982
584 1277.4880318389598
585 1279.3943144383868
586 1279.3943144383873
587 1279.6494688795135
588 1280.9100715253057
589 1280.910071525306
590 1282.6316436464008
591 1282.631643646401
592 1283.2059612838127
593 1283.2415530693595
594 1283.669584127067
595 1283.669584127067
596 1284.916208566405
597 1284.916208566405
598 1285.3576360388322
599 1286.266860423746
600 1286.2668604237463
601 1286.478874868407
602 1289.732190233862
603 1289.921482102238
604 1289.921482102238
605 1292.1006368482347
606 1292.1006368482347
607 1295.6374944409342
608 1306.5417627964537
609 1309.8423097707307
610 1309.8423097707307
611 1309.91124832008
612 1319.607055420309
613 1319.6070554203095
614 1321.5154330069372
615 1321.5154330069377
616 1322.9155295842736
617 1324.5918494233615
618 1324.9704393629622
619 1324.9704393629622
620 1325.9741442683119
621 1325.9741442683119
622 1329.1369875151468
623 1330.429873746001
624 1330.4298737460017
625 1330.691337492324
626 1332.4835382009007
627 1332.5173064546532
628 1332.5173064546536
629 1333.3055982565975
630 1335.4044761791495
631 1339.1663025039618
632 1339.1663025039618
633 1339.4507294657371
634 1340.9426540362197
635 1340.9426540362203
636 1341.6445919103205
637 1343.7592630745544
638 1343.7592630745548
639 1344.876982263801
640 1349.1545463309164
641 1351.07554331032
642 1351.0755433103202
643 1352.4058895240291
644 1352.4058895240291
645 1356.837942233382
Time spent at write_info 0.9959 s
Successfully written nv216_ddh_phonon.xml.
Change the values of <phonon mode>, <asr>, <deltax>, <deltae>,
<ngrid> and <cell> if needed.
Total time required (s): 1.9154143333435059
The code lists all frequencies in cm-1. We see first three acoustic modes are very small (almost 0). Let us check whether the nv216_ddh_phonon.xml file is created or not within the PhonoPy_DynMatrix/ directory.
[3]:
%%bash
ls -lrth PhonoPy_DynMatrix/
total 24M
-rw-rw-r-- 1 arpan arpan 9.2M Dec 21 14:23 FORCE_CONSTANTS
-rw-rw-r-- 1 arpan arpan 8.3K Dec 21 14:23 NV216-DDH-opt.xyz
-rw-rw-r-- 1 arpan arpan 14M Dec 21 14:28 nv216_ddh_phonon.xml
2.2 Generating Stochastically Displaced Configurations
We see that the file is created which will be our input for pyepfd for a frozen-phonon calculation or stochastic calculation. Here we will describe stochastic thermal line sampling method. In PyEPFD this method is named as One Shot Random Sampling (OSR) and we will use the antithetic pairs (AP). Therefore our stochastic algorithm is ‘OSRAP’.
Below we have a small python script osrap.py that takes two command line arguments, first one the phonon file that we just created from FORCE_CONSTANTS and a temperature (in K). Lets look at this python script. You can also use the script in jupyter note book.
[ ]:
# %load osrap.py
import sys,os
from pyepfd.coord_util import *
from pyepfd.pyepfd_io import *
#usage python3 osrap.py phonon_info_file T
phonon_file=sys.argv[1] # First command line argument; phonon file (.xml)
T=sys.argv[2] # Second command line argument; Temperature
# Below we store important phonon file data into sd_inp class.
# From this class we will have access to atoms, coordinates, dynamical matrix,
# acoustic sum rule (asr) etc.
sd_inp = read_pyepfd_info(file_path=phonon_file)
#Below we define an xyz output class using the full file path
#Note the file name would be 300K.xyz if the chosen T = 300
sd_xyz = xyz(file_path=str(T)+'K.xyz',io='w',atoms=sd_inp.atoms)
# Now we will move using stochastic displacement (algo = 'osrap')
sdmoves = ionic_mover( atoms = sd_inp.atoms, \
opt_coord = sd_inp.coord, \
mode = 'SD', \
algo ='osrap', \
ngrid = 500,\
dynmat = sd_inp.ref_dynmatrix, \
asr = sd_inp.asr, \
temperature = float(T) )
# After performing all displacements, we save them into the xyz trajectory file.
for j in range(sdmoves.disp_coord.shape[1]):
sd_xyz.write(cell = sd_inp.cell,coord = sdmoves.disp_coord[:,j])
In the previous file the most important step is the “sdmoves =” lines, specially choosing the mode, algo and ngrid. ngrid is the number of independent sample. Since we are using antithetic pairs for each sample, therefore for ngrid = 500, we will have 1000 configurations. First two configurations would be the antithetic pair of sample 1, then configurations 3 and 4 would be the antithetic pairs of sample 2, and so on.
Now we can run this using a simple slurm job submission script. This code is MPI-parallelized.
Below is the job submission script osrap.sbatch. However if you do not want to submit it in your cluster, you can just run the above python script on your desktop (for small sytems < 100 atoms), it should be fast, for a system of ~200 atoms and ngrid = 500, it takes about 5-10 minutes on a single core (without MPI parallelization).
[ ]:
!/bin/bash
#SBATCH --job-name=ddh_d3h
#SBATCH --time=24:00:00 # MAXIMUM WALLTIME
#SBATCH --partition=gagalli-ivyb # GALLI GROUP PRIVATE PARTITION
#SBATCH --nodes=8
#SBATCH --ntasks-per-node=20
##SBATCH --qos=gagalli-debug
############################## Load distribution ######################################################
export I_MPI_PMI_LIBRARY=/software/slurm-current-el7-x86_64/lib/libpmi.so
export I_MPI_FABRICS=shm:dapl
export OMP_NUM_THREADS=1
export PYTHONUNBUFFERED=TRUE
module load python/anaconda-2020.02 intel/18.0.5 intelmpi/2018.4.274+intel-18.0.5
source activate /home/arpank/Environments/pyepfd
export PATH=/home/arpank/Scripts/pyepfd/utils/qbox_utils:$PATH
phonon_file=PhonoPy_DynMatrix/nv216_ddh_phonon.xml
for T in {0..300..100}; do
## For each OSRAP run we ask 2 nodes and 2*20 cores as each ivyb node has 2 cores
srun --exclusive -N 2 -n 40 python osrap.py ${phonon_file} ${T}> ${T}K.log 2>&1 &
sleep 15
done
wait
This would create four xyz files for each temperature and four log files.
[5]:
%%bash
ls *.xyz *.log
0K.log
0K.xyz
100K.log
100K.xyz
200K.log
200K.xyz
300K.log
300K.xyz
Let us check logfile 300K.log which prints times spent in different parts and if any error occured.
[6]:
%%bash
cat 300K.log
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2024-12-21 14:39:15.178708
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
Process-id0: Time spent on read_pyepfd_info class: 0.527268648147583 s.
Time spent on xyz class: 0.002596616744995117 s.
ionic_mover._stoch_disp(osrap,500) is running with 40 MPI processes.
Process-id = 0: 8% of SD iterations completed.
Process-id = 0: 17% of SD iterations completed.
Process-id = 0: 25% of SD iterations completed.
Process-id = 0: 33% of SD iterations completed.
Process-id = 0: 42% of SD iterations completed.
Process-id = 0: 50% of SD iterations completed.
Process-id = 0: 58% of SD iterations completed.
Process-id = 0: 67% of SD iterations completed.
Process-id = 0: 75% of SD iterations completed.
Process-id = 0: 83% of SD iterations completed.
Process-id = 0: 92% of SD iterations completed.
Process-id = 0: 100% of SD iterations completed.
Time spent on ionic_mover class: 17.660663843154907 s.
300K.xyz is a large file; so we will not open it. But lets check first few lines.
[7]:
%%bash
head 300K.xyz
215
# CELL(abcABC): 10.65000 10.65000 10.65000 90.00000 90.00000 90.00000 PyEPFD-Step: 0 positions{angstrom} cell{angstrom}
C 9.3702659 9.3389617 9.3707882
C 10.176013 10.267503 10.186679
C 9.327185 0.43918429 0.4584821
C 10.258335 1.3662454 1.2811848
C 0.45955569 9.3493541 0.43033944
C 1.3857927 10.250111 1.2786121
C 0.42728748 0.4818096 9.3697941
C 1.3322117 1.3436494 10.192557
As we know first line of xyz file tells us the number of atoms, in this case we have 215 atoms. Second line lists the cell parameters; here we follow a format i-PI code uses. From third line onwards it has cartesian coordinates of 215 atoms. So for each configurations we need 217 lines. Now let us check the total number of lines 300K.xyz has.
[8]:
%%bash
wc -l 300K.xyz
217000 300K.xyz
#So we have 217000/217 = 1000 configurations as it should for ngrid = 500 with antithetic pairs. We save these XYZ files. All 500 sample may not be needed for convergence. Therefore first we would prepare 100 samples, i.e first 100 * 2 = 200 configurations for DFT calculations using quantum espresso. If we get converged results then we would not need to perform DFT calculations for other samples (configurations 201 to 1000).
2.3 Converting XYZ files into a Quantum Espresso Input Deck
To convert that we have a small script xyz2qe.py. Let us have a look at that one.
[ ]:
# %load xyz2qe.py
from pyepfd.coord_util import qe
xyz2qe = qe(io='w', #'w' means writing qe input file mode
path='300K.xyz', # Displaced configuration file
pw_opt_path='TDDFT/pwopts.in', # A file containing QE input options except the coordinates
frames=(1,200), # A tuple containing Range (inclusive) of frames, starting from 1.
pw_path='PWINPs/') # The path where al pw input files would be saved.
This requires 4 input, the xyz file containing displaced coordinates, 300K.xyz, that we just created, the range of frames which in this case is 1 to 200. Both 1 and 200 are included. We need a file that contains all PW options, i.e., functional, cutoffs etc. except the cell parameters and the coordinates. Let us have a look at this file.
[9]:
%%bash
cat TDDFT/pwopts.in
&CONTROL
calculation = 'scf'
restart_mode = 'from_scratch'
verbosity = 'high'
wf_collect = .true.
pseudo_dir='/home/arpank/pseudopotentials/'
tprnfor = .true.
/
&SYSTEM
ibrav = 0
ecutwfc = 60
tot_charge = -1
nspin = 2
tot_magnetization = 2.0
nbnd = 480
nat = 215
ntyp = 2
input_dft = 'PBE0'
exx_fraction = 0.18
ecutfock = 120
/
&ELECTRONS
conv_thr = 1D-8
diago_full_acc=.true.
mixing_beta = 0.3
/
&IONS
ion_dynamics='bfgs'
ion_positions='default'
/
K_POINTS gamma
ATOMIC_SPECIES
C 12.01099968 C_ONCV_PBE-1.0.upf
N 14.00699997 N_ONCV_PBE-1.0.upf
Now we also require a directory where we will save pw input files for each frames. Lets make this directory PWINPs
[10]:
%%bash
mkdir PWINPs
Now let’s convert xyz files into pw input files
[11]:
# %load xyz2qe.py
from pyepfd.coord_util import qe
xyz2qe = qe(io='w', #'w' means writing qe input file mode
path='300K.xyz', # Displaced configuration file
pw_opt_path='TDDFT/pwopts.in', # A file containing QE input options except the coordinates
frames=(1,200), # A tuple containing Range (inclusive) of frames, starting from 1.
pw_path='PWINPs/') # The path where al pw input files would be saved.
del xyz2qe # This line is needed while using it within a jupyter-note book
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2024-12-21 14:54:46.700998
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
Time spent on xyz class: 0.1908116340637207 s.
Time spent on qe class: 0.2591373920440674 s.
Now if we check the PWINPs directory we should be able to see 200 pw input files.
[12]:
%%bash
ls -lrth PWINPs/
total 3.2M
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw8.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw7.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw6.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw5.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw4.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw3.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw2.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw1.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw9.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw21.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw20.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw19.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw18.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw17.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw16.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw15.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw14.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw13.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw12.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw11.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw10.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw32.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw31.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw30.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw29.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw28.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw27.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw26.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw25.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw24.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw23.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw22.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw44.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw43.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw42.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw41.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw40.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw39.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw38.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw37.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw36.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw35.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw34.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw33.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw57.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw56.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw55.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw54.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw53.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw52.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw51.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw50.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw49.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw48.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw47.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw46.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw45.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw69.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw68.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw67.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw66.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw65.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw64.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw63.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw62.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw61.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw60.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw59.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw58.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw82.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw81.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw80.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw79.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw78.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw77.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw76.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw75.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw74.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw73.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw72.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw71.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw70.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw94.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw93.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw92.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw91.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw90.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw89.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw88.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw87.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw86.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw85.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw84.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw83.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw99.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw98.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw97.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw96.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw95.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw107.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw106.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw105.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw104.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw103.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw102.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw101.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw100.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw120.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw119.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw118.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw117.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw116.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw115.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw114.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw113.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw112.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw111.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw110.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw109.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw108.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw132.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw131.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw130.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw129.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw128.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw127.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw126.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw125.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw124.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw123.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw122.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw121.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw144.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw143.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw142.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw141.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw140.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw139.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw138.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw137.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw136.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw135.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw134.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw133.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw157.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw156.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw155.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw154.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw153.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw152.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw151.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw150.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw149.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw148.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw147.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw146.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw145.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw170.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw169.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw168.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw167.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw166.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw165.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw164.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw163.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw162.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw161.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw160.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw159.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw158.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw182.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw181.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw180.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw179.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw178.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw177.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw176.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw175.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw174.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw173.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw172.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw171.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw195.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw194.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw193.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw192.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw191.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw190.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw189.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw188.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw187.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw186.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw185.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw184.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw183.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw199.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw198.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw197.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw196.in
-rw-rw-r-- 1 arpan arpan 13K Dec 21 14:54 pw200.in
Let us have a look at the first 50 lines of such an input file: pw1.in which is for configuration 1.
[13]:
%%bash
head -n 50 PWINPs/pw1.in
&CONTROL
calculation = 'scf'
restart_mode = 'from_scratch'
verbosity = 'high'
wf_collect = .true.
pseudo_dir='/home/arpank/pseudopotentials/'
tprnfor = .true.
/
&SYSTEM
ibrav = 0
ecutwfc = 60
tot_charge = -1
nspin = 2
tot_magnetization = 2.0
nbnd = 480
nat = 215
ntyp = 2
input_dft = 'PBE0'
exx_fraction = 0.18
ecutfock = 120
/
&ELECTRONS
conv_thr = 1D-8
diago_full_acc=.true.
mixing_beta = 0.3
/
&IONS
ion_dynamics='bfgs'
ion_positions='default'
/
K_POINTS gamma
ATOMIC_SPECIES
C 12.01099968 C_ONCV_PBE-1.0.upf
N 14.00699997 N_ONCV_PBE-1.0.upf
CELL_PARAMETERS angstrom
10.65000000 0.00000000 0.00000000
0.00000000 10.65000000 0.00000000
0.00000000 0.00000000 10.65000000
ATOMIC_POSITIONS angstrom
C 9.37026590 9.33896170 9.37078820
C 10.17601300 10.26750300 10.18667900
C 9.32718500 0.43918429 0.45848210
C 10.25833500 1.36624540 1.28118480
C 0.45955569 9.34935410 0.43033944
C 1.38579270 10.25011100 1.27861210
C 0.42728748 0.48180960 9.36979410
C 1.33221170 1.34364940 10.19255700
C 9.39483370 9.38684920 2.24298110
C 10.12114300 10.21578300 3.16942570
C 9.28312160 0.46638730 3.97092770
Now we would submit QE calculations for these input files using a submission script like below.
[ ]:
# %load TDDFT/job.sh
#!/bin/bash
#SBATCH --time=12:00:00
#SBATCH --partition=gagalli-csl
#SBATCH --qos=gagalli-small
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=40
##SBATCH --array=0-7
export nframe=25
module load intel/19.1.1
module load intelmpi/2019.up7+intel-19.1.1
module load mkl/2020.up1
module load python/cpython-3.8.5
#export I_MPI_PMI_LIBRARY=/software/slurm-current-$DISTARCH/lib/libpmi.so
export LD_LIBRARY_PATH=/software/python-3.8.5-el7-x86_64/lib:$LD_LIBRARY_PATH
export OMP_NUM_THREADS=1
QEDIR=/project2/gagalli/jinyu/WEST-Develop/WEST-5.4.0-TDDFT/QEdir12/bin
start_time=$(date +%s)
let begin=151 #${nframe}*${SLURM_ARRAY_TASK_ID}+1
let end=200 #${nframe}*${SLURM_ARRAY_TASK_ID}+${nframe}
workdir=`pwd`
for (( i=$begin; i<=$end; i++ )); do
#mkdir frame-${i}
cd frame-${i}
#mpirun -n 80 ${QEDIR}/pw.x -nb 2 < ${workdir}/PWINPs/pw${i}.in > pw.out
#mpirun -n 40 ${QEDIR}/wbse_init.x -ni 2 -i ${workdir}/wbse_init.in > wbse_init.out
#mpirun -n 80 ${QEDIR}/wbse.x -nb 2 -i ${workdir}/wbse_singlet.in > wbse_singlet.out
#mv pwscf.wbse.save singlet.wbse.save
#mpirun -n 80 ${QEDIR}/wbse.x -nb 2 -i ${workdir}/wbse_triplet.in > wbse_triplet.out
#mv pwscf.wbse.save triplet.wbse.save
mpirun -n 40 ${QEDIR}/westpp.x -i ${workdir}/westpp_singlet.in > westpp_singlet.out
mpirun -n 40 ${QEDIR}/westpp.x -i ${workdir}/westpp_triplet.in > westpp_triplet.out
cd ${workdir}
done
end_time=$(date +%s)
total_time=$((end_time - start_time))
echo "Total execution time: $total_time seconds"
Please remember to load all required modules needed to run quantum espresso on the cluster you are using. Modify those line of the submission scripts accordingly. Using the commented out lines we would perform the further TDDFT calculations using WEST code, afterwards. You have to adapt it according to your need.
2.4 Phonon Renormalizations Single-Particle Levels
Once the calculations are finished we are in a position to analyze them. First let us calculate the phonon-renormalizations of the single-particle defect levels. For that purpose we have to extract the band energies we are interested. For NV-center in diamond we are interesed in bands 430-432, which are the defect levels. Lets extract them using the python script showed below.
[16]:
# %load extract_data.py
import numpy as np
from pyepfd.coord_util import qe
simulation = qe(path = 'TDDFT', frames = (1,200),
bands = (430,432), spin =2)
cell = simulation.cell
atoms = simulation.atoms
coords = simulation.coords
forces = simulation.forces
etotals = simulation.etotals
# Saving an .npz file for future reference
np.savez('NV216_DDH_300K.npz',
cell = cell, atoms = atoms, coords = coords,
forces = forces, etotals = etotals)
del simulation #Needed only within jupyter note book so that file writing is finished.
Time spent on qe class: 0.6937637329101562 s.
Using the above code (extract_data.py), we will extract eigenvalues from both alpha and beta spin channels for the orbitals with indices 430, 431 and 432. All the DFT calculations are saved within he directory DFT. Below is the script using qe_eigval function to perform the following task.
This created two .dat files let us have a look at them.
[17]:
%%bash
ls -lrth *.dat
-rw-rw-r-- 1 arpan arpan 9.9K Dec 21 14:59 is_1_orb-430-432.dat
-rw-rw-r-- 1 arpan arpan 9.9K Dec 21 14:59 is_0_orb-430-432.dat
[18]:
%%bash
head is_0_orb-430-432.dat
#Orbital-430 #Orbital-431 #Orbital-432
13.25855121 14.36966986 14.49676322
13.23186241 14.44880321 14.55506014
13.39142652 14.33276704 14.46798350
13.41372030 14.29533759 14.70180016
13.32512575 14.28971873 14.41910814
13.38637122 14.24825322 14.47878783
13.43827011 14.33585399 14.58418207
13.33799617 14.22218409 14.52373868
13.38736084 14.31245655 14.35497381
Now we will use this files to compute the average and standard deviation of mean using the following script. For that purpose we also need the eigenvalues for the optimized geometry (in the below script we are calling it as a static calculation) as well, which we copied in the directory DFT. The file names are static_is_0.dat and static_is_1.dat, respectively. They contains band energies from 1st to 450 th bands and are extracted using the above script except defining bands = (1,450) and instead of using
open(‘is_’+ …), we used open(‘static_is_’+ …)
[1]:
# %load conv_renorm_levels.py
from pyepfd.epfd import mc_convergence
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import numpy as np
plt.rcParams['font.size'] = 12
orb = [i for i in range(430,433)]
osrap_path = './'
static_path='TDDFT/'
e_scale = 2.0 # Error bar scale, 2 times of standard error
colors = get_cmap('tab10').colors
labels = [[r'$a_1$',r'$e_\mathrm{low}$', r'$e_\mathrm{high}$'],
[r'$\overline{a}_1$',
r'$\overline{e}_\mathrm{low}$',
r'$\overline{e}_\mathrm{high}$']]
label_loc = [[(90,15),(90,-70),(90,150)],
[(90, -50), (90, -135), (90,100)]]
norb = len(orb)
osrap = []; static = [] #; osap_energies=[]
for ispin in range(2):
osrap.append(mc_convergence(\
file_path=osrap_path +'is_' + str(ispin)+ '_orb-' + str(orb[0]) + '-' + str(orb[-1]) + '.dat',
algo='osrap'))
static_orbitals = np.genfromtxt(static_path + '/static_is_' + str(ispin) + '.dat')
static.append([static_orbitals[orb[i]-1] for i in range(norb)])
static = np.array(static)
ndata = len(osrap[0].avg[:,0])
mc = np.zeros((2,norb,ndata),np.float64)
emc = np.zeros((2,norb,ndata),np.float64)
#osap = (osap_energies - static)*1000
for ispin in range(2):
for iorb in range(norb):
mc[ispin,iorb,:] = (osrap[ispin].avg[:,iorb] - static[ispin,iorb])*1000
emc[ispin, iorb,:] = (e_scale*osrap[ispin].sd_mean[:,iorb])*1000
#Printing converged values
outfile = open('conv_renorm_levels.out','w+')
outfile.write("#Converged Mean renorm. values after %d steps\n"%ndata)
outfile.write("#-------------------------------------------------------------------\n")
outfile.write("#Ispin Band-index Mean_Renorm(meV) %3.1f*SD_Mean(meV)\n"%e_scale)
outfile.write("#-------------------------------------------------------------------\n")
for ispin in range(2):
for iorb in range(norb):
outfile.write(" %2d %5d %12.4f %12.4f\n"
%(ispin, orb[iorb], mc[ispin,iorb,-1], emc[ispin, iorb,-1]))
outfile.close()
fig, axs = plt.subplots(1,2)
#fig = []; axs = []
for ispin in range(2):
#tmpfig, tmpaxs = plt.subplots(num='ispin='+str(ispin))
#fig.append(tmpfig); axs.append(tmpaxs)
axs[ispin].set_xlabel("Number of samples")
axs[ispin].set_ylabel("Renormalization (meV)")
for iorb in range(norb):
axs[ispin].plot([i+1 for i in range(ndata)], mc[ispin,iorb,:],
color=colors[iorb], label = labels[ispin][iorb])
axs[ispin].fill_between([i+1 for i in range(ndata)],\
mc[ispin,iorb,:] + emc[ispin,iorb,:], mc[ispin,iorb,:] - emc[ispin,iorb,:], alpha = 0.3, color=colors[iorb])
#axs[ispin].axhline(osap[ispin,iorb],color=colors[iorb],linestyle='dashed')
axs[ispin].text(label_loc[ispin][iorb][0], label_loc[ispin][iorb][1],
labels[ispin][iorb], color=colors[iorb],
ha='center', va='center')
#axs[ispin].legend()
#for ispin in range(2):
# fig[ispin].savefig(fig_name+'_isp'+str(ispin)+'.png',
# dpi=600,bbox_inches = 'tight', pad_inches = 0.1)
axs[0].set_title(r'$\alpha$-spin channel')
axs[1].set_title(r'$\beta$-spin channel')
fig.set_size_inches(6.6,3.3)
plt.subplots_adjust(wspace=0.5)
plt.show()
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2024-12-21 17:39:56.619589
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
The above scripts also prints the renormalization values into the follwing files:
[2]:
%%bash
cat conv_renorm_levels.out
#Converged Mean renorm. values after 100 steps
#-------------------------------------------------------------------
#Ispin Band-index Mean_Renorm(meV) 2.0*SD_Mean(meV)
#-------------------------------------------------------------------
0 430 52.4587 12.1295
0 431 -121.2145 13.7529
0 432 114.1844 15.1193
1 430 -1.5932 18.9691
1 431 -173.6870 14.8185
1 432 19.0860 14.9976
2.5 Phonon Renormalizations of Many-Body States
Once we have the DFT calculations and obtained the renormalization of the single particle levels; we can start doing TDDFT calculations for each frames using the WEST code. This can be performed adjusting the previous job submission script by commenting out the appropriate lines.
We need the following file the wbse_init calculation only for hybrid TDDFT (not required for GGA functionals). Check the WEST documentation: https://west-code.org/doc/West/latest/, for more details.
[ ]:
# %load TDDFT/wbse_init.in
input_west:
qe_prefix: pwscf
west_prefix: pwscf
outdir: ./
wbse_init_control:
wbse_init_calculation: S
solver: TDDFT
localization: W
overlap_thr: 0.001
We need the following input file: wbse_singlet.in for spin-flip TDDFT afterwards.
[ ]:
# %load TDDFT/wbse_singlet.in
input_west:
qe_prefix: pwscf
west_prefix: pwscf
outdir: ./
wbse_init_control:
wbse_init_calculation: S
solver: TDDFT
localization: W
overlap_thr: 0.001
wbse_control:
wbse_calculation: D
scissor_ope: 0.0
n_liouville_eigen: 4
n_liouville_times: 10
n_liouville_maxiter: 1000
n_liouville_read_from_file: 0
trev_liouville: 0.000000001
trev_liouville_rel: 0.0000001
wbse_epsinfty: 1.0
spin_excitation: S
l_preconditioning: True
l_pre_shift: True
l_spin_flip: True
l_spin_flip_kernel: True
l_spin_flip_alda0: False
l_print_spin_flip_kernel: False
spin_flip_cut1: 500
spin_flip_cut2: 0.0001
l_reduce_io: True
l_minimize_exx_if_active: False
n_exx_lowrank: 480
While we need the following input file: wbse_triplet.in for spin-conserved TDDFT
[ ]:
# %load TDDFT/wbse_triplet.in
input_west:
qe_prefix: pwscf
west_prefix: pwscf
outdir: ./
wbse_init_control:
wbse_init_calculation: S
solver: TDDFT
localization: W
overlap_thr: 0.001
wbse_control:
wbse_calculation: D
scissor_ope: 0.0
n_liouville_eigen: 2
n_liouville_times: 20
n_liouville_maxiter: 1000
n_liouville_read_from_file: 0
trev_liouville: 0.000000001
trev_liouville_rel: 0.0000001
wbse_epsinfty: 1.0
spin_excitation: S
l_preconditioning: True
l_pre_shift: False
l_reduce_io: True
l_minimize_exx_if_active: False
n_exx_lowrank: 480
Once all the calculations for singlet and triplet states are finished, we can start analyzing the data. First step is to extract the renormalizations of the singlet and triplet states for each frame. we can do that using the following script.
[3]:
# %load tddft_renorm.py
import sys,os
import numpy as np
from pyepfd.west_util import tddft_energy
from pyepfd.constants import *
ev2mev = 1000
# Energies of the optimized GS geometry
singlet_0 = np.array([
0.11141222885296939E-1,
0.60939419639036742E-1,
0.60945316008269398E-1,
0.15486512004137221E+0
]) / ev2unit['Ry']
triplet_0 = np.array([
0.17576134841668695E+0,
0.1757629043174222E+0
]) / ev2unit['Ry']
# Obtaining tddft energies of each frame
singlet = tddft_energy( path = 'TDDFT',
west_prefix = 'singlet',
frames = (1,200) )
triplet = tddft_energy( path = 'TDDFT',
west_prefix = 'triplet',
frames = (1,200) )
## singlet renormalizations
# The first value is related to ms=0 state of the triplet state;
# therefore we are subtracting this value as it should be 0.
singlet_e_low = (singlet[:,1] - singlet[:,0] - singlet_0[1] + singlet_0[0]) * ev2mev
singlet_e_high = (singlet[:,2] - singlet[:,0] - singlet_0[2] + singlet_0[0]) * ev2mev
singlet_a = (singlet[:,3] - singlet[:,0] - singlet_0[3] + singlet_0[0]) * ev2mev
## triplet renormalizations for each frames
triplet_e_low = (triplet[:,0] - triplet_0[0]) * ev2mev
triplet_e_high = (triplet[:,1] - triplet_0[1]) * ev2mev
## writing data in files
singlet_out = open('singlet_renorm.dat','w+')
triplet_out = open('triplet_renorm.dat','w+')
singlet_out.write('# 1E(low) 1E(High) 1A1 all in meV\n')
triplet_out.write('# 3E(low) 3E(High) all in meV\n')
for frame in range(len(singlet_e_low)):
singlet_out.write(' %10.2f %10.2f %10.2f\n'\
%(singlet_e_low[frame],singlet_e_high[frame],singlet_a[frame]))
triplet_out.write(' %10.2f %10.2f \n'\
%(triplet_e_low[frame],triplet_e_high[frame]))
singlet_out.close()
triplet_out.close()
singlet: TDDFT energies for frames: (1, 200) extracted in 0.009213685989379883 s.
triplet: TDDFT energies for frames: (1, 200) extracted in 0.009042501449584961 s.
Lets have a look at the newly created singlet_renorm.dat and triplet_renorm.dat files containg phonon renormalizations of the many-body states of interest.
[4]:
%%bash
echo "================================================="
head -n 10 singlet_renorm.dat
echo "================================================="
head -n 10 triplet_renorm.dat
=================================================
# 1E(low) 1E(High) 1A1 all in meV
-113.15 193.16 74.11
-109.57 98.72 -18.52
-124.97 45.09 -36.30
-80.53 7.55 57.09
-50.54 7.41 -63.71
-128.61 41.71 -76.50
-19.62 12.40 -6.08
-102.35 -2.66 -78.31
-49.75 17.09 -38.91
=================================================
# 3E(low) 3E(High) all in meV
165.05 379.24
54.24 262.40
-260.47 -131.53
-180.81 113.70
-142.70 -20.44
-249.09 -87.45
32.83 257.54
-210.48 -3.69
-147.42 -9.11
Now we will use this files to compute the average and standard deviation of mean using the following script.
[9]:
# %load conv_renorm_states.py
from pyepfd.epfd import mc_convergence
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import numpy as np
plt.rcParams['font.size'] = 12
fig_name='excited_DDH_nv216_osrap_300K'
colors = get_cmap('tab10').colors
e_scale = 2.0
label_loc = [[(90,-108),(90,120),(90,0)],
[(90, -100), (90, 80)]]
class tddft_osrap_conv:
def __init__(self,path):
self.renorm = \
mc_convergence(file_path=path, algo='osrap')
self.ndata = len(self.renorm.avg[:,0])
self.nstates = np.shape(self.renorm.avg)[1]
if self.nstates == 3:
self.states = ['1E(low)','1E(high)','1A1']
self.labels = [r'$^{1}E_\mathrm{low}$',
r'$^{1}E_\mathrm{high}$',
r'$^{1}A_{1}$']
elif self.nstates == 2:
self.states = ['3E(low)','3E(high)']
self.labels = [r'$^{3}E_\mathrm{low}$',
r'$^{3}E_\mathrm{high}$']
singlet = tddft_osrap_conv("singlet_renorm.dat")
triplet = tddft_osrap_conv("triplet_renorm.dat")
all_states = [singlet, triplet]
#Printing converged values
outfile = open('conv_renorm_many_body.out','w+')
outfile.write("#Converged Mean renorm. values after %d steps\n"%singlet.ndata)
outfile.write("#-------------------------------------------------------------------\n")
outfile.write("#State Mean_Renorm(meV) %3.1f*SD_Mean(meV)\n"%e_scale)
outfile.write("#-------------------------------------------------------------------\n")
for state in all_states:
for istate in range(state.nstates):
outfile.write(" %8s %12.4f %12.4f\n"
%(state.states[istate], state.renorm.avg[-1,istate],
state.renorm.sd_mean[-1,istate]*e_scale))
outfile.close()
fig, axs = plt.subplots(1,2,figsize=(6.6,3.3))
for k in range(len(all_states)):
axs[k].set_xlabel("Number of samples")
axs[k].set_ylabel("Renormalization (meV)")
for j in range(all_states[k].nstates):
axs[k].plot([i+1 for i in range(all_states[k].ndata)],
all_states[k].renorm.avg[:,j], label = all_states[k].labels[j],
color = colors[j])
axs[k].fill_between([i+1 for i in range(all_states[k].ndata)],
all_states[k].renorm.avg[:,j] + all_states[k].renorm.sd_mean[:,j]*e_scale,
all_states[k].renorm.avg[:,j] - all_states[k].renorm.sd_mean[:,j]*e_scale,
alpha = 0.3, color = colors[j])
#axs[k].set_xlim(0,200)
#axs[k].set_ylim(-150,100)
#axs[k].legend()
axs[k].text(label_loc[k][j][0], label_loc[k][j][1],
all_states[k].labels[j], color=colors[j],
ha='center', va='center')
axs[0].set_title(r'Singlet excited states')
axs[1].set_title(r'Triplet excited states')
fig.set_size_inches(6.6,3.3)
plt.subplots_adjust(wspace=0.5)
#fig.savefig(fig_name+'.png',dpi=600,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
The above scripts also prints the renormalization values into the follwing file:
[8]:
%%bash
cat conv_renorm_many_body.out
#Converged Mean renorm. values after 100 steps
#-------------------------------------------------------------------
#State Mean_Renorm(meV) 2.0*SD_Mean(meV)
#-------------------------------------------------------------------
1E(low) -142.7186 9.2154
1E(high) 83.3722 10.9285
1A1 -27.9857 7.6659
3E(low) -154.0781 22.5723
3E(high) 18.7320 22.3219
The analysis scripts are available at the directory: pyepfd_path/tutorials/2_NV. Due to the large sizes, TDDFT/frame-n directories, PhonoPy_DynMatrix directory and xyz files files are removed from the PyEPFD distribution.
These results are published at the following paper:
A.Kundu, G. Galli, “Quantum Vibronic Effects on the Excitation Energies of the Nitrogen-Vacancy Center in Diamond”, J. Phys. Chem. Lett., 15, 802 (2024). https://pubs.acs.org/doi/10.1021/acs.jpclett.3c03269
Arxiv Link: https://arxiv.org/abs/2401.06745