MICCoM School 2022 Ex#3 : More systems (input generation with ASE)¶
You can try different systems, here are some examples:
Benzene molecule
H\(_2\)O
SiC solid
We will use the ASE (Atomic Simulation Environment) Python package to build the input files.
3.1 : Benzene molecule¶
[1]:
# download the pseudopotential files
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.2.upf
[2]:
# build the input for DFT
from ase import Atom, Atoms
from ase.cell import Cell
from ase.io import write
cell = Cell([[10,0,0],[0,10,0],[0,0,10]]) # Angstroms
positions = []
positions.append(Atom('C',[1.2115,0.6995,0.0000])) # Angstroms
positions.append(Atom('C',[1.2115,-0.6995,0.0000]))
positions.append(Atom('C',[0.0000,1.3990,0.0000]))
positions.append(Atom('C',[0.0000,-1.3990,0.0000]))
positions.append(Atom('C',[-1.2115,-0.6995,0.0000]))
positions.append(Atom('C',[-1.2115,0.6995,0.0000]))
positions.append(Atom('H',[0.0000,2.5000,0.0000]))
positions.append(Atom('H',[2.1651,1.2500,0.0000]))
positions.append(Atom('H',[2.1651,-1.2500,0.0000]))
positions.append(Atom('H',[0.0000,-2.5000,0.0000]))
positions.append(Atom('H',[-2.1651,-1.2500,0.0000]))
positions.append(Atom('H',[-2.1651,1.2500,0.0000]))
atoms = Atoms(positions,cell=cell)
pseudopotentials = {}
pseudopotentials['C'] = 'C_ONCV_PBE-1.2.upf'
pseudopotentials['H'] = 'H_ONCV_PBE-1.2.upf'
inp_data = {}
inp_data['calculation'] = 'scf'
inp_data['restart_mode'] = 'from_scratch'
inp_data['prefix'] = 'benzene'
inp_data['outdir'] = './'
inp_data['pseudo_dir'] = './'
inp_data['verbosity'] = 'high'
inp_data['wf_collect'] = True
inp_data['ecutwfc'] = 50
inp_data['nbnd'] = 30
inp_data['assume_isolated'] = 'mp'
inp_data['diago_full_acc'] = True
write('benzene.pwi',atoms,pseudopotentials=pseudopotentials,input_data=inp_data)
[3]:
!cat benzene.pwi
&CONTROL
calculation = 'scf'
verbosity = 'high'
restart_mode = 'from_scratch'
wf_collect = .true.
outdir = './'
prefix = 'benzene'
pseudo_dir = './'
/
&SYSTEM
nbnd = 30
ecutwfc = 50
assume_isolated = 'mp'
ntyp = 2
nat = 12
ibrav = 0
/
&ELECTRONS
diago_full_acc = .true.
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
C 12.011 C_ONCV_PBE-1.2.upf
H 1.008 H_ONCV_PBE-1.2.upf
K_POINTS gamma
CELL_PARAMETERS angstrom
10.00000000000000 0.00000000000000 0.00000000000000
0.00000000000000 10.00000000000000 0.00000000000000
0.00000000000000 0.00000000000000 10.00000000000000
ATOMIC_POSITIONS angstrom
C 1.2115000000 0.6995000000 0.0000000000
C 1.2115000000 -0.6995000000 0.0000000000
C 0.0000000000 1.3990000000 0.0000000000
C 0.0000000000 -1.3990000000 0.0000000000
C -1.2115000000 -0.6995000000 0.0000000000
C -1.2115000000 0.6995000000 0.0000000000
H 0.0000000000 2.5000000000 0.0000000000
H 2.1651000000 1.2500000000 0.0000000000
H 2.1651000000 -1.2500000000 0.0000000000
H 0.0000000000 -2.5000000000 0.0000000000
H -2.1651000000 -1.2500000000 0.0000000000
H -2.1651000000 1.2500000000 0.0000000000
[4]:
# build the input for MBPT
import yaml
west = {}
west['input_west'] = {}
west['input_west']['qe_prefix'] = 'benzene'
west['input_west']['west_prefix'] = 'benzene'
west['input_west']['outdir'] = './'
west['wstat_control'] = {}
west['wstat_control']['wstat_calculation'] = 'S'
west['wstat_control']['n_pdep_eigen'] = 50
west['wstat_control']['n_steps_write_restart'] = 0
west['wfreq_control'] = {}
west['wfreq_control']['wfreq_calculation'] = 'XWGQ'
west['wfreq_control']['n_pdep_eigen_to_use'] = 50
west['wfreq_control']['qp_bandrange'] = [13,18]
west['wfreq_control']['o_restart_time'] = -1.0
with open('west.in', 'w') as file :
yaml.dump(west, file, sort_keys=False)
[5]:
!cat west.in
input_west:
qe_prefix: benzene
west_prefix: benzene
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
n_steps_write_restart: 0
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 50
qp_bandrange:
- 13
- 18
o_restart_time: -1.0
Run pw.x, wstat.x, and wfreq.x. These calculations should finish in 5 minutes.
[6]:
!mpirun -n 8 pw.x -i benzene.pwi > benzene.pwo
!mpirun -n 64 wstat.x -ni 8 -i west.in > wstat.out
!mpirun -n 64 wfreq.x -ni 8 -i west.in > wfreq.out
3.2 : H2O¶
[7]:
# download the pseudopotential files
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/O_ONCV_PBE-1.2.upf
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.2.upf
!wget -N -q http://west-code.org/doc/training/water/H2O.cif
[8]:
# build the input for DFT, starting from a CIF file
from ase.io import read, write
atoms = read('H2O.cif')
pseudopotentials = {}
pseudopotentials['O'] = 'O_ONCV_PBE-1.2.upf'
pseudopotentials['H'] = 'H_ONCV_PBE-1.2.upf'
inp_data = {}
inp_data['calculation'] = 'scf'
inp_data['restart_mode'] = 'from_scratch'
inp_data['prefix'] = 'h2o'
inp_data['outdir'] = './'
inp_data['pseudo_dir'] = './'
inp_data['verbosity'] = 'high'
inp_data['wf_collect'] = True
inp_data['ecutwfc'] = 50
inp_data['nbnd'] = 80
inp_data['diago_full_acc'] = True
write('h2o.pwi',atoms,pseudopotentials=pseudopotentials,input_data=inp_data)
[9]:
!cat h2o.pwi
&CONTROL
calculation = 'scf'
verbosity = 'high'
restart_mode = 'from_scratch'
wf_collect = .true.
outdir = './'
prefix = 'h2o'
pseudo_dir = './'
/
&SYSTEM
nbnd = 80
ecutwfc = 50
ntyp = 2
nat = 36
ibrav = 0
/
&ELECTRONS
diago_full_acc = .true.
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
H 1.008 H_ONCV_PBE-1.2.upf
O 15.999 O_ONCV_PBE-1.2.upf
K_POINTS gamma
CELL_PARAMETERS angstrom
7.60356630000000 0.00000000000000 0.00000000000000
-3.80178315000000 6.58488157515925 0.00000000000000
0.00000000000000 0.00000000000000 7.14296200000000
ATOMIC_POSITIONS angstrom
H 2.5154346141 0.0000000000 1.4030777397
H -1.2577173070 2.1784302773 1.4030777397
H 2.5440658430 4.4064512978 1.4030777397
H 5.0881316859 0.0000000000 4.9745587397
H -2.5440658430 4.4064512978 4.9745587397
H 1.2577173070 2.1784302773 4.9745587397
H 3.4394504051 0.0000000000 0.0760439735
H -1.7197252026 2.9786514259 0.0760439735
H 2.0820579474 3.6062301493 0.0760439735
H 4.1641158949 0.0000000000 3.6475249735
H -2.0820579474 3.6062301493 3.6475249735
H 1.7197252026 2.9786514259 3.6475249735
H 1.7103728160 5.7785299019 6.9918669248
H -2.0577569496 5.1768429328 6.9918669248
H 0.3473841335 2.2143903156 6.9918669248
H 5.5121559660 0.8063516733 6.9918669248
H 1.7440262004 1.4080386423 6.9918669248
H 4.1491672835 4.3704912595 6.9918669248
H 2.0914103340 0.8063516733 3.4203859248
H 5.8595400996 1.4080386423 3.4203859248
H 3.4543990165 4.3704912595 3.4203859248
H -1.7103728160 5.7785299019 3.4203859248
H 2.0577569496 5.1768429328 3.4203859248
H -0.3473841335 2.2143903156 3.4203859248
O 2.4836441034 0.0000000000 0.3975629790
O -1.2418220517 2.1508988875 0.3975629790
O 2.5599610983 4.4339826877 0.3975629790
O 5.1199221966 0.0000000000 3.9690439790
O -2.5599610983 4.4339826877 3.9690439790
O 1.2418220517 2.1508988875 3.9690439790
O 5.0229539156 0.0000000000 6.6399189011
O -2.5114769578 4.3500056930 6.6399189011
O 1.2903061922 2.2348758822 6.6399189011
O 2.5806123844 0.0000000000 3.0684379011
O -1.2903061922 2.2348758822 3.0684379011
O 2.5114769578 4.3500056930 3.0684379011
[10]:
# build the input for MBPT
import yaml
west = {}
west['input_west'] = {}
west['input_west']['qe_prefix'] = 'h2o'
west['input_west']['west_prefix'] = 'h2o'
west['input_west']['outdir'] = './'
west['wstat_control'] = {}
west['wstat_control']['wstat_calculation'] = 'S'
west['wstat_control']['n_pdep_eigen'] = 200
west['wstat_control']['n_steps_write_restart'] = 0
west['wfreq_control'] = {}
west['wfreq_control']['wfreq_calculation'] = 'XWGQ'
west['wfreq_control']['macropol_calculation'] = 'C'
# Note that for a solid we must set macropol_calculation to 'C', the default 'N' is good for molecules
west['wfreq_control']['n_pdep_eigen_to_use'] = 200
west['wfreq_control']['qp_bandrange'] = [20,30]
west['wfreq_control']['o_restart_time'] = -1.0
with open('west.in', 'w') as file :
yaml.dump(west, file, sort_keys=False)
[11]:
!cat west.in
input_west:
qe_prefix: h2o
west_prefix: h2o
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 200
n_steps_write_restart: 0
wfreq_control:
wfreq_calculation: XWGQ
macropol_calculation: C
n_pdep_eigen_to_use: 200
qp_bandrange:
- 20
- 30
o_restart_time: -1.0
Run pw.x, wstat.x, and wfreq.x. These calculations should finish in 10 minutes.
[12]:
!mpirun -n 8 pw.x -i h2o.pwi > h2o.pwo
!mpirun -n 64 wstat.x -ni 8 -i west.in > wstat.out
!mpirun -n 64 wfreq.x -ni 8 -i west.in > wfreq.out
3.3 : SiC solid¶
[13]:
# download the pseudopotential files
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/Si_ONCV_PBE-1.2.upf
[14]:
# build the input for DFT
from ase.build import bulk
from ase.io import write
atoms = bulk('SiC', 'zincblende', a=4.36, cubic=True) # Angstrom
pseudopotentials = {}
pseudopotentials['C'] = 'C_ONCV_PBE-1.2.upf'
pseudopotentials['Si'] = 'Si_ONCV_PBE-1.2.upf'
inp_data = {}
inp_data['calculation'] = 'scf'
inp_data['restart_mode'] = 'from_scratch'
inp_data['prefix'] = 'sic'
inp_data['outdir'] = './'
inp_data['pseudo_dir'] = './'
inp_data['verbosity'] = 'high'
inp_data['wf_collect'] = True
inp_data['ecutwfc'] = 30
write('sic_scf.pwi',atoms,pseudopotentials=pseudopotentials,input_data=inp_data,kpts=[8,8,8])
inp_data['calculation'] = 'nscf'
inp_data['nosym'] = True
inp_data['noinv'] = True
inp_data['nbnd'] = 20
inp_data['diago_full_acc'] = True
# Note that for systems with k-points we run a nscf calculation making sure that symmetries are turned off
write('sic_nscf.pwi',atoms,pseudopotentials=pseudopotentials,input_data=inp_data,kpts=[2,2,2])
[15]:
!cat sic_scf.pwi
!cat sic_nscf.pwi
&CONTROL
calculation = 'scf'
verbosity = 'high'
restart_mode = 'from_scratch'
wf_collect = .true.
outdir = './'
prefix = 'sic'
pseudo_dir = './'
/
&SYSTEM
ecutwfc = 30
ntyp = 2
nat = 8
ibrav = 0
/
&ELECTRONS
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
Si 28.085 Si_ONCV_PBE-1.2.upf
C 12.011 C_ONCV_PBE-1.2.upf
K_POINTS automatic
8 8 8 0 0 0
CELL_PARAMETERS angstrom
4.36000000000000 0.00000000000000 0.00000000000000
0.00000000000000 4.36000000000000 0.00000000000000
0.00000000000000 0.00000000000000 4.36000000000000
ATOMIC_POSITIONS angstrom
Si 0.0000000000 0.0000000000 0.0000000000
C 1.0900000000 1.0900000000 1.0900000000
Si 0.0000000000 2.1800000000 2.1800000000
C 1.0900000000 3.2700000000 3.2700000000
Si 2.1800000000 0.0000000000 2.1800000000
C 3.2700000000 1.0900000000 3.2700000000
Si 2.1800000000 2.1800000000 0.0000000000
C 3.2700000000 3.2700000000 1.0900000000
&CONTROL
calculation = 'nscf'
verbosity = 'high'
restart_mode = 'from_scratch'
wf_collect = .true.
outdir = './'
prefix = 'sic'
pseudo_dir = './'
/
&SYSTEM
nbnd = 20
ecutwfc = 30
nosym = .true.
noinv = .true.
ntyp = 2
nat = 8
ibrav = 0
/
&ELECTRONS
diago_full_acc = .true.
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
Si 28.085 Si_ONCV_PBE-1.2.upf
C 12.011 C_ONCV_PBE-1.2.upf
K_POINTS automatic
2 2 2 0 0 0
CELL_PARAMETERS angstrom
4.36000000000000 0.00000000000000 0.00000000000000
0.00000000000000 4.36000000000000 0.00000000000000
0.00000000000000 0.00000000000000 4.36000000000000
ATOMIC_POSITIONS angstrom
Si 0.0000000000 0.0000000000 0.0000000000
C 1.0900000000 1.0900000000 1.0900000000
Si 0.0000000000 2.1800000000 2.1800000000
C 1.0900000000 3.2700000000 3.2700000000
Si 2.1800000000 0.0000000000 2.1800000000
C 3.2700000000 1.0900000000 3.2700000000
Si 2.1800000000 2.1800000000 0.0000000000
C 3.2700000000 3.2700000000 1.0900000000
[16]:
# build the input for MBPT
import yaml
west = {}
west['input_west'] = {}
west['input_west']['qe_prefix'] = 'sic'
west['input_west']['west_prefix'] = 'sic'
west['input_west']['outdir'] = './'
west['wstat_control'] = {}
west['wstat_control']['wstat_calculation'] = 'S'
west['wstat_control']['n_pdep_eigen'] = 48
west['wstat_control']['n_steps_write_restart'] = 0
west['wfreq_control'] = {}
west['wfreq_control']['wfreq_calculation'] = 'XWGQ'
west['wfreq_control']['macropol_calculation'] = 'C'
west['wfreq_control']['n_pdep_eigen_to_use'] = 48
west['wfreq_control']['qp_bandrange'] = [15,18]
west['wfreq_control']['o_restart_time'] = -1.0
with open('west.in', 'w') as file :
yaml.dump(west, file, sort_keys=False)
[17]:
!cat west.in
input_west:
qe_prefix: sic
west_prefix: sic
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 48
n_steps_write_restart: 0
wfreq_control:
wfreq_calculation: XWGQ
macropol_calculation: C
n_pdep_eigen_to_use: 48
qp_bandrange:
- 15
- 18
o_restart_time: -1.0
Run pw.x, wstat.x, and wfreq.x. These calculations should finish in 15 minutes.
[18]:
!mpirun -n 16 pw.x -i sic_scf.pwi > sic_scf.pwo
!mpirun -n 16 pw.x -i sic_nscf.pwi > sic_nscf.pwo
!mpirun -n 64 wstat.x -ni 8 -i west.in > wstat.out
!mpirun -n 64 wfreq.x -ni 8 -i west.in > wfreq.out