Skip to content

Commit c8b6d25

Browse files
Merge branch 'master' into Diegom_sobol
2 parents 909a6bb + 1c12b51 commit c8b6d25

32 files changed

Lines changed: 7755 additions & 640 deletions
Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
import cellconstructor as CC, cellconstructor.Phonons, cellconstructor.ForceTensor
2+
import ase, ase.dft.kpoints
3+
4+
import numpy as np
5+
import matplotlib.pyplot as plt
6+
import matplotlib.cm as cm
7+
import matplotlib.colors as colors
8+
9+
import sys, os
10+
11+
NQIRR = 13
12+
#CMAP = "Spectral_r"
13+
PATH = "GXWXKGL"
14+
N_POINTS = 1000
15+
16+
SPECIAL_POINTS = {"G": [0,0,0],
17+
"X": [0, .5, .5],
18+
"L": [.5, .5, .5],
19+
"W": [.25, .75, .5],
20+
"K": [3/8., 3/4., 3/8.]}
21+
22+
# Load the harmonic and sscha phonons
23+
harmonic_dyn = CC.Phonons.Phonons('harmonic_dyn', NQIRR)
24+
sscha_dyn = CC.Phonons.Phonons('sscha_T300_dyn', NQIRR)
25+
26+
# Get the band path
27+
qpath, data = CC.Methods.get_bandpath(harmonic_dyn.structure.unit_cell,
28+
PATH,
29+
SPECIAL_POINTS,
30+
N_POINTS)
31+
xaxis, xticks, xlabels = data # Info to plot correclty the x axis
32+
33+
# Get the phonon dispersion along the path
34+
harmonic_dispersion = CC.ForceTensor.get_phonons_in_qpath(harmonic_dyn, qpath)
35+
sscha_dispersion = CC.ForceTensor.get_phonons_in_qpath(sscha_dyn, qpath)
36+
37+
nmodes = harmonic_dyn.structure.N_atoms * 3
38+
39+
# Plot the two dispersions
40+
plt.figure(dpi = 150)
41+
ax = plt.gca()
42+
43+
for i in range(nmodes):
44+
lbl=None
45+
lblsscha = None
46+
if i == 0:
47+
lbl = 'Harmonic'
48+
lblsscha = 'SSCHA'
49+
50+
ax.plot(xaxis, harmonic_dispersion[:,i], color = 'k', ls = 'dashed', label = lbl)
51+
ax.plot(xaxis, sscha_dispersion[:,i], color = 'r', label = lblsscha)
52+
53+
# Plot vertical lines for each high symmetry points
54+
for x in xticks:
55+
ax.axvline(x, 0, 1, color = "k", lw = 0.4)
56+
ax.axhline(0, 0, 1, color = 'k', ls = ':', lw = 0.4)
57+
58+
ax.legend()
59+
60+
# Set the x labels to the high symmetry points
61+
ax.set_xticks(xticks)
62+
ax.set_xticklabels(xlabels)
63+
64+
ax.set_xlabel("Q path")
65+
ax.set_ylabel("Phonons [cm-1]")
66+
67+
plt.tight_layout()
68+
plt.savefig("dispersion.png")
69+
plt.show()
70+
71+
72+

Examples/ThermodynamicsOfGold/plot_thermal_expansion.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@
5252

5353
# Prepare the figure and plot the V(T) from the sscha data
5454
plt.figure(dpi = 150)
55-
plt.scatter(temperatures, volumes, label = "SSCHA data")
55+
plt.scatter(temperatures, volumes, label = "SSCHA data", color = 'r')
5656

5757
# Fit the data with a quadratic curve
5858
def parabola(x, a, b, c):
@@ -72,12 +72,12 @@ def diff_parab(x, a, b, c):
7272

7373
# Plot the fit
7474
_t_ = np.linspace(np.min(temperatures), np.max(temperatures), 1000)
75-
plt.plot(_t_, parabola(_t_, *popt), label = "Fit")
75+
plt.plot(_t_, parabola(_t_, *popt), label = "Fit", color = 'k', zorder = 0)
7676

7777
# Adjust the plot adding labels, legend, and saving in eps
7878
plt.xlabel("Temperature [K]")
7979
plt.ylabel(r"Volume [$\AA^3$]")
8080
plt.legend()
8181
plt.tight_layout()
82-
plt.savefig("thermal_expansion.eps")
82+
plt.savefig("thermal_expansion.png")
8383
plt.show()

Examples/sscha_and_dft/Au_ONCV_PBE-1.0.oncvpsp.upf

Lines changed: 5049 additions & 0 deletions
Large diffs are not rendered by default.
Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
# Import the sscha code
2+
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities
3+
4+
# Import the cellconstructor library to manage phonons
5+
import cellconstructor as CC, cellconstructor.Phonons
6+
import cellconstructor.Structure, cellconstructor.calculators
7+
8+
# Import the DFT calculator
9+
import cellconstructor.calculators
10+
11+
# Import numerical and general pourpouse libraries
12+
import numpy as np, matplotlib.pyplot as plt
13+
import sys, os
14+
15+
16+
# Initialize the DFT (Quantum Espresso) calculator for gold
17+
# The input data is a dictionary that encodes the pw.x input file namelist
18+
input_data = {
19+
'control' : {
20+
# Avoid writing wavefunctions on the disk
21+
'disk_io' : 'None',
22+
# Where to find the pseudopotential
23+
'pseudo_dir' : '.'
24+
},
25+
'system' : {
26+
# Specify the basis set cutoffs
27+
'ecutwfc' : 45, # Cutoff for wavefunction
28+
'ecutrho' : 45*4, # Cutoff for the density
29+
# Information about smearing (it is a metal)
30+
'occupations' : 'smearing',
31+
'smearing' : 'mv',
32+
'degauss' : 0.03
33+
},
34+
'electrons' : {
35+
'conv_thr' : 1e-8
36+
}
37+
}
38+
39+
# the pseudopotential for each chemical element
40+
# In this case just Gold
41+
pseudopotentials = {'Au' : 'Au_ONCV_PBE-1.0.oncvpsp.upf'}
42+
43+
# the kpoints mesh and the offset
44+
kpts = (1,1,1)
45+
koffset = (1,1,1)
46+
47+
command = 'mpirun -np 4 pw.x -npool 1 -i PREFIX.pwi > PREFIX.pwo'
48+
49+
# Prepare the quantum espresso calculator
50+
calculator = CC.calculators.Espresso(input_data,
51+
pseudopotentials,
52+
command = command,
53+
kpts = kpts,
54+
koffset = koffset)
55+
56+
57+
58+
TEMPERATURE = 300
59+
N_CONFIGS = 50
60+
MAX_ITERATIONS = 20
61+
START_DYN = 'start_dyn'
62+
NQIRR = 13
63+
64+
# Let us load the starting dynamical matrix
65+
gold_dyn = CC.Phonons.Phonons(START_DYN, NQIRR)
66+
67+
# Initialize the random ionic ensemble
68+
ensemble = sscha.Ensemble.Ensemble(gold_dyn, TEMPERATURE)
69+
70+
# Initialize the free energy minimizer
71+
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
72+
minim.set_minimization_step(0.01)
73+
74+
# Initialize the NVT simulation
75+
relax = sscha.Relax.SSCHA(minim, calculator, N_configs = N_CONFIGS,
76+
max_pop = MAX_ITERATIONS)
77+
78+
# Define the I/O operations
79+
# To save info about the free energy minimization after each step
80+
ioinfo = sscha.Utilities.IOInfo()
81+
ioinfo.SetupSaving("minim_info")
82+
relax.setup_custom_functions(custom_function_post = ioinfo.CFP_SaveAll)
83+
84+
85+
# Run the NVT simulation (save the stress to compute the pressure)
86+
relax.relax(get_stress = True)
87+
88+
# If instead you want to run a NPT simulation, use
89+
# The target pressure is given in GPa.
90+
#relax.vc_relax(target_press = 0)
91+
92+
# You can also run a mixed simulation (NVT) but with variable lattice parameters
93+
#relax.vc_relax(fix_volume = True)
94+
95+
# Now we can save the final dynamical matrix
96+
# And print in stdout the info about the minimization
97+
relax.minim.finalize()
98+
relax.minim.dyn.save_qe("sscha_T{}_dyn".format(TEMPERATURE))
99+
100+
101+
102+
103+
104+
105+
106+
107+
108+

Examples/sscha_and_dft/start_dyn1

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
Dynamical matrix file
2+
File generated with the CellConstructor by Lorenzo Monacelli
3+
1 1 0 1.8897259890000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
4+
Basis vectors
5+
2.9495460300000000 0.0000000000000000 0.0000000000000000
6+
1.4747730150000002 2.5543817916115379 0.0000000000000000
7+
1.4747730150000002 0.8514605972038461 2.4082942487839478
8+
1 'Au ' 179524.0513667391205672
9+
1 1 0.0000000000000000 0.0000000000000000 0.0000000000000000
10+
11+
Dynamical Matrix in cartesian axes
12+
13+
q = ( 0.000000000000 0.000000000000 0.000000000000 )
14+
15+
1 1
16+
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
17+
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
18+
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
19+
20+
Diagonalizing the dynamical matrix
21+
22+
q = ( 0.000000000000 0.000000000000 0.000000000000 )
23+
24+
***************************************************************************
25+
freq ( 1) = 0.00000000 [THz] = 0.00000000 [cm-1]
26+
( 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 )
27+
freq ( 2) = 0.00000000 [THz] = 0.00000000 [cm-1]
28+
( 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 )
29+
freq ( 3) = 0.00000000 [THz] = 0.00000000 [cm-1]
30+
( 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 )
31+
***************************************************************************

Examples/sscha_and_dft/start_dyn10

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
Dynamical matrix file
2+
File generated with the CellConstructor by Lorenzo Monacelli
3+
1 1 0 1.8897259890000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
4+
Basis vectors
5+
2.9495460300000000 0.0000000000000000 0.0000000000000000
6+
1.4747730150000002 2.5543817916115379 0.0000000000000000
7+
1.4747730150000002 0.8514605972038461 2.4082942487839478
8+
1 'Au ' 179524.0513667391205672
9+
1 1 0.0000000000000000 0.0000000000000000 0.0000000000000000
10+
11+
Dynamical Matrix in cartesian axes
12+
13+
q = ( 0.169517612173 0.097871039020 0.069205275373 )
14+
15+
1 1
16+
0.0751622289125264 0.0000000000000000 0.0281657539230274 0.0000000000000000 0.0209028729909581 0.0000000000000000
17+
0.0281657539230274 0.0000000000000000 0.0426392176937492 0.0000000000000000 0.0120682793481662 0.0000000000000000
18+
0.0209028729909581 0.0000000000000000 0.0120682793481662 0.0000000000000000 0.0353854466021494 0.0000000000000000
19+
20+
Dynamical Matrix in cartesian axes
21+
22+
q = ( 0.169517612173 -0.097871039020 -0.069205275373 )
23+
24+
1 1
25+
0.0751622289125264 0.0000000000000000 -0.0281657539230274 0.0000000000000000 -0.0209028729909581 0.0000000000000000
26+
-0.0281657539230274 0.0000000000000000 0.0426392176937492 0.0000000000000000 0.0120682793481662 0.0000000000000000
27+
-0.0209028729909581 0.0000000000000000 0.0120682793481662 0.0000000000000000 0.0353854466021494 0.0000000000000000
28+
29+
Dynamical Matrix in cartesian axes
30+
31+
q = ( 0.000000000000 -0.195742078041 0.069205275373 )
32+
33+
1 1
34+
0.0263777120843606 0.0000000000000000 -0.0000000000000000 0.0000000000000000 -0.0000000000000000 0.0000000000000000
35+
-0.0000000000000000 0.0000000000000000 0.0914237345219151 0.0000000000000000 -0.0241365586963325 0.0000000000000000
36+
-0.0000000000000000 0.0000000000000000 -0.0241365586963325 0.0000000000000000 0.0353854466021494 0.0000000000000000
37+
38+
Diagonalizing the dynamical matrix
39+
40+
q = ( 0.169517612173 0.097871039020 0.069205275373 )
41+
42+
***************************************************************************
43+
freq ( 1) = 1.26105061 [THz] = 42.06411449 [cm-1]
44+
( -0.500000 0.000000 0.866025 0.000000 0.000000 0.000000 )
45+
freq ( 2) = 1.26213016 [THz] = 42.10012438 [cm-1]
46+
( -0.301466 0.000000 -0.174051 0.000000 0.937456 0.000000 )
47+
freq ( 3) = 2.46009088 [THz] = 82.05978700 [cm-1]
48+
( -0.811861 0.000000 -0.468728 0.000000 -0.348103 0.000000 )
49+
***************************************************************************

Examples/sscha_and_dft/start_dyn11

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
Dynamical matrix file
2+
File generated with the CellConstructor by Lorenzo Monacelli
3+
1 1 0 1.8897259890000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
4+
Basis vectors
5+
2.9495460300000000 0.0000000000000000 0.0000000000000000
6+
1.4747730150000002 2.5543817916115379 0.0000000000000000
7+
1.4747730150000002 0.8514605972038461 2.4082942487839478
8+
1 'Au ' 179524.0513667391205672
9+
1 1 0.0000000000000000 0.0000000000000000 0.0000000000000000
10+
11+
Dynamical Matrix in cartesian axes
12+
13+
q = ( -0.169517612173 -0.097871039020 0.138410550746 )
14+
15+
1 1
16+
0.0817675664305969 0.0000000000000000 0.0142701466013224 0.0000000000000000 -0.0183567296612535 0.0000000000000000
17+
0.0142701466013224 0.0000000000000000 0.0652898204672991 0.0000000000000000 -0.0105982628113659 0.0000000000000000
18+
-0.0183567296612535 0.0000000000000000 -0.0105982628113659 0.0000000000000000 0.0755111590679396 0.0000000000000000
19+
20+
Dynamical Matrix in cartesian axes
21+
22+
q = ( -0.169517612173 0.097871039020 -0.138410550746 )
23+
24+
1 1
25+
0.0817675664305969 0.0000000000000000 -0.0142701466013224 0.0000000000000000 0.0183567296612535 0.0000000000000000
26+
-0.0142701466013224 0.0000000000000000 0.0652898204672991 0.0000000000000000 -0.0105982628113659 0.0000000000000000
27+
0.0183567296612535 0.0000000000000000 -0.0105982628113659 0.0000000000000000 0.0755111590679395 0.0000000000000000
28+
29+
Dynamical Matrix in cartesian axes
30+
31+
q = ( 0.000000000000 0.195742078041 0.138410550746 )
32+
33+
1 1
34+
0.0570509474856502 0.0000000000000000 -0.0000000000000000 0.0000000000000000 -0.0000000000000000 0.0000000000000000
35+
-0.0000000000000000 0.0000000000000000 0.0900064394122459 0.0000000000000000 0.0211965256227318 0.0000000000000000
36+
-0.0000000000000000 0.0000000000000000 0.0211965256227317 0.0000000000000000 0.0755111590679396 0.0000000000000000
37+
38+
Diagonalizing the dynamical matrix
39+
40+
q = ( -0.169517612173 -0.097871039020 0.138410550746 )
41+
42+
***************************************************************************
43+
freq ( 1) = 1.85457920 [THz] = 61.86209446 [cm-1]
44+
( 0.500000 0.000000 -0.866025 0.000000 0.000000 0.000000 )
45+
freq ( 2) = 1.90756503 [THz] = 63.62951127 [cm-1]
46+
( 0.503661 0.000000 0.290789 0.000000 0.813491 0.000000 )
47+
freq ( 3) = 2.51790636 [THz] = 83.98830356 [cm-1]
48+
( -0.704504 0.000000 -0.406746 0.000000 0.581577 0.000000 )
49+
***************************************************************************

Examples/sscha_and_dft/start_dyn12

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
Dynamical matrix file
2+
File generated with the CellConstructor by Lorenzo Monacelli
3+
1 1 0 1.8897259890000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
4+
Basis vectors
5+
2.9495460300000000 0.0000000000000000 0.0000000000000000
6+
1.4747730150000002 2.5543817916115379 0.0000000000000000
7+
1.4747730150000002 0.8514605972038461 2.4082942487839478
8+
1 'Au ' 179524.0513667391205672
9+
1 1 0.0000000000000000 0.0000000000000000 0.0000000000000000
10+
11+
Dynamical Matrix in cartesian axes
12+
13+
q = ( 0.000000000000 0.000000000000 0.103807913060 )
14+
15+
1 1
16+
0.0137058449993487 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
17+
-0.0000000000000000 0.0000000000000000 0.0137058449993487 0.0000000000000000 -0.0000000000000000 0.0000000000000000
18+
0.0000000000000000 0.0000000000000000 -0.0000000000000000 0.0000000000000000 0.0888425361262248 0.0000000000000000
19+
20+
Dynamical Matrix in cartesian axes
21+
22+
q = ( 0.000000000000 0.000000000000 -0.103807913060 )
23+
24+
1 1
25+
0.0137058449993487 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
26+
-0.0000000000000000 0.0000000000000000 0.0137058449993487 0.0000000000000000 -0.0000000000000000 0.0000000000000000
27+
0.0000000000000000 0.0000000000000000 -0.0000000000000000 0.0000000000000000 0.0888425361262248 0.0000000000000000
28+
29+
Diagonalizing the dynamical matrix
30+
31+
q = ( 0.000000000000 0.000000000000 0.103807913060 )
32+
33+
***************************************************************************
34+
freq ( 1) = 0.90900624 [THz] = 30.32117990 [cm-1]
35+
( 0.000000 0.000000 1.000000 0.000000 -0.000000 0.000000 )
36+
freq ( 2) = 0.90900624 [THz] = 30.32117990 [cm-1]
37+
( 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 )
38+
freq ( 3) = 2.31432543 [THz] = 77.19757585 [cm-1]
39+
( 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 )
40+
***************************************************************************

0 commit comments

Comments
 (0)