cp2k: implement full cell spec parsing

* previously we only allowed for orthorombic cells, with the help of the
  new cp2k-input-tools we can now parse any cell spec
* switching to cp2k-input-tools allows to parse almost any CP2K input
* switch to using Angstrom (CP2K's default for distances)
This commit is contained in:
Tiziano Müller 2019-09-16 09:32:39 +02:00
parent c14d77d032
commit 88c3b12129
8 changed files with 251 additions and 380 deletions

View File

@ -1,20 +0,0 @@
This is an example of CP2K interface.
To create supercells with displacements:
% phonopy --cp2k -c Si.inp -d --dim="2 2 2"
A perfect 2x2x2 supercell (supercell.in) and one 2x2x2 supercells
(supercell-xxx.in) of the conventional unit cell written in Si.in are
created. In addition, disp.yaml file is created. After force
calculation with the crystal structure in supercell-001.in, it is
needed to create FORCE_SETS file by
% phonopy --cp2k -f Si-supercell-001-forces-1_0.xyz
Here .out files are supposed to contain the forces on atoms calculated
by CP2K and filenames can be chosen freely. The disp.yaml file has
to be put in the current directory. Now you can run phonon
calculation, e.g.,
% phonopy --cp2k -c Si.inp -p --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/2 1/2 0 0 0 1/2 0 1/2"

26
example/Si-CP2K/README.md Normal file
View File

@ -0,0 +1,26 @@
# Example for the CP2K Phonopy interface using bulk silicon
To create supercells with displacements:
```console
$ phonopy --cp2k -c Si.inp -d --dim="2 2 2"
```
A perfect 2x2x2 supercell (`Si-supercell-000.inp`) and one 2x2x2 supercells
(`supercell-001.inp`) of the conventional unit cell written in `Si.inp` are
created. In addition, a `phonopy_disp.yaml` file is created. After the force
calculation with the crystal structure in `supercell-001.inp`, it is
needed to create `FORCE_SETS` file by running:
```console
$ phonopy --cp2k -f Si-supercell-001-forces-1_0.xyz
```
Here the `.xyz` files are supposed to contain the forces on atoms calculated
by CP2K.
To plot the phonon band structure:
```console
$ phonopy --cp2k -c Si.inp -p --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/2 1/2 0 0 0 1/2 0 1/2"
```

View File

@ -2,68 +2,68 @@
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 Si -0.00273314 -0.00000000 -0.00000000
2 1 Si 0.00000278 -0.00000000 -0.00000000
3 1 Si 0.00001016 -0.00000000 -0.00000000
4 1 Si 0.00001060 -0.00000000 0.00000000
5 1 Si 0.00001016 -0.00000000 -0.00000000
6 1 Si 0.00001060 -0.00000000 0.00000000
7 1 Si -0.00006163 -0.00000000 -0.00000000
8 1 Si 0.00000535 -0.00000000 -0.00000000
9 1 Si -0.00010297 0.00002371 0.00002370
10 1 Si -0.00000227 0.00000123 0.00000123
11 1 Si -0.00010290 0.00002332 -0.00002332
12 1 Si -0.00000223 0.00000123 -0.00000123
13 1 Si -0.00010290 -0.00002332 0.00002332
14 1 Si -0.00000223 -0.00000123 0.00000123
15 1 Si -0.00010297 -0.00002370 -0.00002370
16 1 Si -0.00000227 -0.00000123 -0.00000123
17 1 Si 0.00003110 -0.00002346 0.00004240
18 1 Si 0.00003137 -0.00002377 -0.00004277
19 1 Si 0.00000152 -0.00000112 -0.00001093
20 1 Si 0.00000148 -0.00000112 0.00001091
21 1 Si 0.00003110 0.00002346 -0.00004240
22 1 Si 0.00003137 0.00002377 0.00004277
23 1 Si 0.00000152 0.00000111 0.00001093
24 1 Si 0.00000148 0.00000112 -0.00001091
25 1 Si 0.00003110 0.00004240 -0.00002346
26 1 Si 0.00003137 -0.00004277 -0.00002377
27 1 Si 0.00003110 -0.00004240 0.00002346
28 1 Si 0.00003137 0.00004277 0.00002377
29 1 Si 0.00000152 -0.00001093 -0.00000112
30 1 Si 0.00000148 0.00001091 -0.00000112
31 1 Si 0.00000152 0.00001093 0.00000112
32 1 Si 0.00000148 -0.00001091 0.00000112
33 1 Si -0.00000476 0.00000901 0.00000901
34 1 Si 0.00004197 0.00000643 0.00000643
35 1 Si 0.00000657 0.00000620 0.00000727
36 1 Si -0.00000782 0.00000918 0.00000787
37 1 Si 0.00000657 0.00000727 0.00000620
38 1 Si -0.00000782 0.00000787 0.00000918
39 1 Si 0.00000100 0.00000797 0.00000797
40 1 Si 0.00069379 0.00042676 0.00042676
41 1 Si 0.00000100 -0.00000797 -0.00000797
42 1 Si 0.00069379 -0.00042676 -0.00042676
43 1 Si 0.00000657 -0.00000727 -0.00000620
44 1 Si -0.00000782 -0.00000787 -0.00000918
45 1 Si 0.00000657 -0.00000620 -0.00000727
46 1 Si -0.00000782 -0.00000918 -0.00000787
47 1 Si -0.00000476 -0.00000901 -0.00000901
48 1 Si 0.00004197 -0.00000643 -0.00000643
49 1 Si -0.00000782 -0.00000780 0.00000929
50 1 Si 0.00000668 -0.00000723 0.00000622
51 1 Si 0.00068699 -0.00041420 0.00041420
52 1 Si 0.00000101 -0.00000787 0.00000787
53 1 Si 0.00004152 -0.00000632 0.00000632
54 1 Si -0.00000474 -0.00000896 0.00000896
55 1 Si -0.00000782 -0.00000929 0.00000780
56 1 Si 0.00000668 -0.00000622 0.00000723
57 1 Si -0.00000782 0.00000929 -0.00000780
58 1 Si 0.00000668 0.00000622 -0.00000723
59 1 Si 0.00004152 0.00000632 -0.00000632
60 1 Si -0.00000474 0.00000896 -0.00000896
61 1 Si 0.00068699 0.00041420 -0.00041420
62 1 Si 0.00000101 0.00000787 -0.00000787
63 1 Si -0.00000782 0.00000780 -0.00000929
64 1 Si 0.00000668 0.00000723 -0.00000622
SUM OF ATOMIC FORCES 0.00000002 -0.00000000 -0.00000000 0.00000002
1 1 Si -0.00258246 -0.00000000 0.00000000
2 1 Si 0.00000263 -0.00000000 -0.00000000
3 1 Si 0.00000960 -0.00000000 0.00000000
4 1 Si 0.00001002 -0.00000000 0.00000000
5 1 Si 0.00000960 -0.00000000 -0.00000000
6 1 Si 0.00001002 -0.00000000 0.00000000
7 1 Si -0.00005823 -0.00000000 -0.00000000
8 1 Si 0.00000505 -0.00000000 -0.00000000
9 1 Si -0.00009729 0.00002239 0.00002239
10 1 Si -0.00000214 0.00000116 0.00000116
11 1 Si -0.00009723 0.00002205 -0.00002205
12 1 Si -0.00000210 0.00000116 -0.00000116
13 1 Si -0.00009723 -0.00002205 0.00002205
14 1 Si -0.00000210 -0.00000116 0.00000116
15 1 Si -0.00009729 -0.00002239 -0.00002239
16 1 Si -0.00000214 -0.00000116 -0.00000116
17 1 Si 0.00002940 -0.00002218 0.00004007
18 1 Si 0.00002963 -0.00002245 -0.00004040
19 1 Si 0.00000143 -0.00000105 -0.00001033
20 1 Si 0.00000140 -0.00000106 0.00001031
21 1 Si 0.00002940 0.00002218 -0.00004007
22 1 Si 0.00002963 0.00002245 0.00004040
23 1 Si 0.00000143 0.00000105 0.00001033
24 1 Si 0.00000140 0.00000106 -0.00001031
25 1 Si 0.00002940 0.00004007 -0.00002218
26 1 Si 0.00002963 -0.00004040 -0.00002245
27 1 Si 0.00002940 -0.00004007 0.00002218
28 1 Si 0.00002963 0.00004041 0.00002245
29 1 Si 0.00000143 -0.00001033 -0.00000105
30 1 Si 0.00000140 0.00001031 -0.00000106
31 1 Si 0.00000143 0.00001033 0.00000105
32 1 Si 0.00000140 -0.00001031 0.00000106
33 1 Si -0.00000450 0.00000851 0.00000851
34 1 Si 0.00003964 0.00000607 0.00000607
35 1 Si 0.00000621 0.00000586 0.00000687
36 1 Si -0.00000739 0.00000868 0.00000743
37 1 Si 0.00000621 0.00000687 0.00000586
38 1 Si -0.00000739 0.00000744 0.00000868
39 1 Si 0.00000095 0.00000753 0.00000753
40 1 Si 0.00065536 0.00040290 0.00040290
41 1 Si 0.00000095 -0.00000753 -0.00000753
42 1 Si 0.00065536 -0.00040290 -0.00040290
43 1 Si 0.00000621 -0.00000687 -0.00000586
44 1 Si -0.00000739 -0.00000744 -0.00000868
45 1 Si 0.00000621 -0.00000586 -0.00000687
46 1 Si -0.00000739 -0.00000868 -0.00000743
47 1 Si -0.00000450 -0.00000851 -0.00000851
48 1 Si 0.00003964 -0.00000607 -0.00000607
49 1 Si -0.00000739 -0.00000737 0.00000878
50 1 Si 0.00000630 -0.00000684 0.00000588
51 1 Si 0.00064929 -0.00039169 0.00039169
52 1 Si 0.00000096 -0.00000744 0.00000744
53 1 Si 0.00003924 -0.00000597 0.00000597
54 1 Si -0.00000448 -0.00000847 0.00000847
55 1 Si -0.00000739 -0.00000877 0.00000737
56 1 Si 0.00000630 -0.00000587 0.00000684
57 1 Si -0.00000739 0.00000877 -0.00000737
58 1 Si 0.00000630 0.00000587 -0.00000684
59 1 Si 0.00003924 0.00000597 -0.00000597
60 1 Si -0.00000448 0.00000847 -0.00000847
61 1 Si 0.00064929 0.00039169 -0.00039169
62 1 Si 0.00000096 0.00000744 -0.00000744
63 1 Si -0.00000739 0.00000737 -0.00000878
64 1 Si 0.00000630 0.00000684 -0.00000588
SUM OF ATOMIC FORCES 0.00000002 0.00000000 0.00000000 0.00000002

View File

@ -1,204 +0,0 @@
natom: 64
displacements:
- atom: 1
direction:
[ 1.0000000000000000, 0.0000000000000000, 0.0000000000000000 ]
displacement:
[ 0.0200000000000000, 0.0000000000000000, 0.0000000000000000 ]
lattice:
- [ 20.659105641913271, 0.000000000000000, 0.000000000000000 ] # a
- [ 0.000000000000000, 20.659105641913271, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 20.659105641913271 ] # c
points:
- symbol: Si # 1
coordinates: [ 0.437500000000000, 0.437500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 2
coordinates: [ 0.937500000000000, 0.437500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 3
coordinates: [ 0.437500000000000, 0.937500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 4
coordinates: [ 0.937500000000000, 0.937500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 5
coordinates: [ 0.437500000000000, 0.437500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 6
coordinates: [ 0.937500000000000, 0.437500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 7
coordinates: [ 0.437500000000000, 0.937500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 8
coordinates: [ 0.937500000000000, 0.937500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 9
coordinates: [ 0.437500000000000, 0.187500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 10
coordinates: [ 0.937500000000000, 0.187500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 11
coordinates: [ 0.437500000000000, 0.687500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 12
coordinates: [ 0.937500000000000, 0.687500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 13
coordinates: [ 0.437500000000000, 0.187500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 14
coordinates: [ 0.937500000000000, 0.187500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 15
coordinates: [ 0.437500000000000, 0.687500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 16
coordinates: [ 0.937500000000000, 0.687500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 17
coordinates: [ 0.187500000000000, 0.437500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 18
coordinates: [ 0.687500000000000, 0.437500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 19
coordinates: [ 0.187500000000000, 0.937500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 20
coordinates: [ 0.687500000000000, 0.937500000000000, 0.187500000000000 ]
mass: 28.085500
- symbol: Si # 21
coordinates: [ 0.187500000000000, 0.437500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 22
coordinates: [ 0.687500000000000, 0.437500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 23
coordinates: [ 0.187500000000000, 0.937500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 24
coordinates: [ 0.687500000000000, 0.937500000000000, 0.687500000000000 ]
mass: 28.085500
- symbol: Si # 25
coordinates: [ 0.187500000000000, 0.187500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 26
coordinates: [ 0.687500000000000, 0.187500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 27
coordinates: [ 0.187500000000000, 0.687500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 28
coordinates: [ 0.687500000000000, 0.687500000000000, 0.437500000000000 ]
mass: 28.085500
- symbol: Si # 29
coordinates: [ 0.187500000000000, 0.187500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 30
coordinates: [ 0.687500000000000, 0.187500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 31
coordinates: [ 0.187500000000000, 0.687500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 32
coordinates: [ 0.687500000000000, 0.687500000000000, 0.937500000000000 ]
mass: 28.085500
- symbol: Si # 33
coordinates: [ 0.062500000000000, 0.062500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 34
coordinates: [ 0.562500000000000, 0.062500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 35
coordinates: [ 0.062500000000000, 0.562500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 36
coordinates: [ 0.562500000000000, 0.562500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 37
coordinates: [ 0.062500000000000, 0.062500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 38
coordinates: [ 0.562500000000000, 0.062500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 39
coordinates: [ 0.062500000000000, 0.562500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 40
coordinates: [ 0.562500000000000, 0.562500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 41
coordinates: [ 0.062500000000000, 0.312500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 42
coordinates: [ 0.562500000000000, 0.312500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 43
coordinates: [ 0.062500000000000, 0.812500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 44
coordinates: [ 0.562500000000000, 0.812500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 45
coordinates: [ 0.062500000000000, 0.312500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 46
coordinates: [ 0.562500000000000, 0.312500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 47
coordinates: [ 0.062500000000000, 0.812500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 48
coordinates: [ 0.562500000000000, 0.812500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 49
coordinates: [ 0.312500000000000, 0.062500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 50
coordinates: [ 0.812500000000000, 0.062500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 51
coordinates: [ 0.312500000000000, 0.562500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 52
coordinates: [ 0.812500000000000, 0.562500000000000, 0.312500000000000 ]
mass: 28.085500
- symbol: Si # 53
coordinates: [ 0.312500000000000, 0.062500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 54
coordinates: [ 0.812500000000000, 0.062500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 55
coordinates: [ 0.312500000000000, 0.562500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 56
coordinates: [ 0.812500000000000, 0.562500000000000, 0.812500000000000 ]
mass: 28.085500
- symbol: Si # 57
coordinates: [ 0.312500000000000, 0.312500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 58
coordinates: [ 0.812500000000000, 0.312500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 59
coordinates: [ 0.312500000000000, 0.812500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 60
coordinates: [ 0.812500000000000, 0.812500000000000, 0.062500000000000 ]
mass: 28.085500
- symbol: Si # 61
coordinates: [ 0.312500000000000, 0.312500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 62
coordinates: [ 0.812500000000000, 0.312500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 63
coordinates: [ 0.312500000000000, 0.812500000000000, 0.562500000000000 ]
mass: 28.085500
- symbol: Si # 64
coordinates: [ 0.812500000000000, 0.812500000000000, 0.562500000000000 ]
mass: 28.085500

View File

@ -99,7 +99,9 @@ def write_supercells_with_displacements(interface_mode,
optional_structure_info[1])
elif interface_mode == 'cp2k':
from phonopy.interface.cp2k import write_supercells_with_displacements
write_supercells_with_displacements(supercell, cells_with_disps)
write_supercells_with_displacements(supercell,
cells_with_disps,
optional_structure_info)
elif interface_mode == 'crystal':
from phonopy.interface.crystal import (
write_supercells_with_displacements)
@ -165,8 +167,8 @@ def read_crystal_structure(filename=None,
return unitcell, (cell_filename, atypes)
elif interface_mode == 'cp2k':
from phonopy.interface.cp2k import read_cp2k
unitcell = read_cp2k(cell_filename)
return unitcell, (cell_filename,)
unitcell, config_tree = read_cp2k(cell_filename)
return unitcell, (cell_filename, config_tree)
elif interface_mode == 'crystal':
from phonopy.interface.crystal import read_crystal
unitcell, conv_numbers = read_crystal(cell_filename)
@ -216,7 +218,7 @@ def get_default_supercell_filename(interface_mode):
elif interface_mode == 'siesta':
return "supercell.fdf"
elif interface_mode == 'cp2k':
return "supercell.inp"
return None # CP2K interface generates filenames based on original project name
elif interface_mode == 'crystal':
return None # supercell.ext can not be parsed by crystal interface.
elif interface_mode == 'dftbp':
@ -228,12 +230,9 @@ def get_default_supercell_filename(interface_mode):
def get_default_displacement_distance(interface_mode):
if interface_mode in ('wien2k', 'abinit', 'elk', 'qe', 'siesta', 'cp2k',
'turbomole'):
if interface_mode in ('wien2k', 'abinit', 'elk', 'qe', 'siesta', 'turbomole'):
displacement_distance = 0.02
elif interface_mode == 'crystal':
displacement_distance = 0.01
else: # default or vasp
else: # default or vasp, crystal, cp2k
displacement_distance = 0.01
return displacement_distance
@ -251,6 +250,7 @@ def get_default_physical_units(interface_mode=None):
CRYSTAL : eV, Angstrom, AMU, eV/Angstroem
DFTB+ : hartree, au, AMU hartree/au
TURBOMOLE : hartree, au, AMU, hartree/au
CP2K : hartree, Angstrom, AMU, hartree/au
"""
@ -302,9 +302,9 @@ def get_default_physical_units(interface_mode=None):
units['length_unit'] = 'au'
elif interface_mode == 'cp2k':
units['factor'] = CP2KToTHz
units['nac_factor'] = Hartree / Bohr # in a.u.
units['distance_to_A'] = Bohr
units['force_constants_unit'] = 'hartree/au^2'
units['nac_factor'] = None # not implemented
units['distance_to_A'] = 1.0
units['force_constants_unit'] = 'hartree/Angstrom.au'
units['length_unit'] = 'Angstrom'
elif interface_mode == 'crystal':
units['factor'] = CrystalToTHz
@ -475,7 +475,8 @@ def get_force_constant_conversion_factor(unit, interface_mode):
'eV/Angstrom.au': 1 / Bohr,
'Ry/au^2': Rydberg / Bohr ** 2,
'mRy/au^2': Rydberg / Bohr ** 2 / 1000,
'hartree/au^2': Hartree / Bohr ** 2}
'hartree/au^2': Hartree / Bohr ** 2,
'hartree/Angstrom.au': Hartree / Bohr }
if default_unit not in factor_to_eVperA2:
msg = "Force constant conversion for %s unit is not implemented."
raise NotImplementedError(msg)

View File

@ -1,5 +1,5 @@
# vim: set fileencoding=utf-8 :
# Copyright (C) 2017 Tiziano Müller
# Copyright (C) 2017-2019 Tiziano Müller
# All rights reserved.
#
# This file is part of phonopy.
@ -33,124 +33,191 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import print_function
import sys
from fractions import Fraction
from os import path
import numpy as np
from phonopy.file_IO import collect_forces
from phonopy.file_IO import iter_collect_forces
from phonopy.interface.vasp import (check_forces, get_drift_forces)
from phonopy.structure.atoms import (PhonopyAtoms, symbol_map)
from phonopy.units import Bohr
def parse_set_of_forces(num_atoms, forces_filenames, verbose=True):
hook = '# Atom Kind Element'
force_sets = []
for i, filename in enumerate(forces_filenames):
if verbose:
sys.stdout.write("%d. " % (i + 1))
with open(filename) as fhandle:
cp2k_forces = collect_forces(fhandle, num_atoms, hook, [3, 4, 5])
if check_forces(cp2k_forces, num_atoms, filename, verbose=verbose):
drift_force = get_drift_forces(cp2k_forces,
filename=filename,
verbose=verbose)
force_sets.append(np.array(cp2k_forces) - drift_force)
else:
return [] # if one file is invalid, the whole thing is broken
forces = iter_collect_forces(filename, num_atoms, '# Atom Kind Element', [3, 4, 5])
if not check_forces(forces, num_atoms, filename, verbose=verbose):
return [] # if one file is invalid, the whole thing is broken
drift_force = get_drift_forces(forces, filename=filename, verbose=verbose)
force_sets.append(np.array(forces) - drift_force)
return force_sets
def read_cp2k(filename):
from cp2k_tools.parser import CP2KInputParser
from cp2k_input_tools.parser import CP2KInputParser
with open(filename) as fhandle:
parser = CP2KInputParser()
cp2k_in = parser.parse(fhandle)
tree = parser.parse(fhandle)
cp2k_cell = cp2k_in['force_eval']['subsys']['cell']
cp2k_coord = cp2k_in['force_eval']['subsys']['coord']
try:
subsys = tree['+force_eval'][0]['+subsys']
cp2k_cell = subsys['+cell']
except IndexError:
raise RuntimeError("could not find a FORCE_EVAL/SUBSYS/CELL section in the given CP2K input file")
if 'abc' in cp2k_cell:
unit = '[angstrom]' # CP2K default unit
unit_cell = np.eye(3)
if isinstance(cp2k_cell['abc'][0], str):
unit = cp2k_cell['abc'][0]
unit_cell *= cp2k_cell['abc'][1:]
if len(tree['+force_eval']) > 1:
raise NotImplementedError("the given CP2K input file contains multiple FORCE_EVAL sections, which is not (yet) supported")
# CP2K can get its cell information in two ways:
# - A, B, C: cell vectors
# - ABC: scaling of cell vectors, ALPHA_BETA_GAMMA: angles between the cell vectors
# We'll parse either of them, but only write A, B, C
if 'a' in cp2k_cell:
# unit vectors given
unit_cell = np.array([
cp2k_cell['a'],
cp2k_cell['b'],
cp2k_cell['c'],
])
elif 'abc' in cp2k_cell:
# length of unit vectors given
if 'alpha_beta_gamma' in cp2k_cell:
# if we also have the angles, construct the cell
alpha, beta, gamma = cp2k_cell.pop('alpha_beta_gamma')
cos_alpha = np.cos(alpha)
cos_beta = np.cos(beta)
cos_gamma = np.cos(gamma)
sin_gamma = np.sin(gamma)
unit_cell = np.array([
[1., 0., 0.],
[cos_gamma, sin_gamma, 0.],
[
cos_beta,
(cos_alpha-cos_gamma*cos_beta)/sin_gamma,
np.sqrt(1.-cos_beta**2-((cos_alpha-cos_gamma*cos_beta)/sin_gamma)**2),
],
])
else:
unit_cell *= cp2k_cell['abc']
unit_cell = np.eye(3)
if unit == '[angstrom]':
unit_cell /= Bohr # phonopy expects the lattice to be in Bohr
else:
raise NotImplementedError("unit scaling for other units than angstrom not yet implemented")
a, b, c = cp2k_cell.pop('abc') # remove them from the tree since we pass it along
unit_cell[0,:] *= a
unit_cell[1,:] *= b
unit_cell[2,:] *= c
if '+cell_ref' in cp2k_cell:
print("WARNING: the &CELL_REF section must be manually adjusted")
cp2k_coord = subsys['+coord']
numbers = []
positions = []
for coordline in cp2k_coord['*']:
# coordinates are a series of strings according to the CP2K schema
fields = coordline.split()
numbers += [symbol_map[fields[0]]]
# positions can also be fractions
positions += [[float(Fraction(f)) for f in fields[1:4]]]
if cp2k_coord.get('scaled', False): # the keyword can be unavailable, true or false
return (PhonopyAtoms(numbers=numbers, cell=unit_cell, scaled_positions=positions), tree)
else:
raise NotImplementedError("unit cell can only be specified via ABC")
# the keyword can be unavailable, true, false or None, with unavailable=false, None=true
scaled_coords = cp2k_coord.get('scaled', False) is not False
if not scaled_coords:
raise NotImplementedError("only scaled coordinates are currently supported")
numbers = [symbol_map[e[0]] for e in cp2k_coord['*']]
positions = [e[1:] for e in cp2k_coord['*']]
return PhonopyAtoms(numbers=numbers, cell=unit_cell, scaled_positions=positions)
def write_cp2k(filename, cell):
with open(filename, 'w') as fhandle:
fhandle.write(get_cp2k_structure(cell))
return (PhonopyAtoms(numbers=numbers, cell=unit_cell, positions=positions), tree)
def write_supercells_with_displacements(supercell,
cells_with_displacements,
optional_structure_info,
pre_filename="supercell",
width=3):
write_cp2k("supercell.inp", supercell)
orig_fname, tree = optional_structure_info
fbase, fext = path.splitext(orig_fname)
pbase = tree["+global"]["project_name"]
supercell_ref_name = "{}-supercell{}".format(fbase, fext)
with open(supercell_ref_name, 'w') as fhandle:
fhandle.write("""\
# Generated by Phonopy, based on {fname}
# Original configuration with the generated supercell for comparison
""".format(fname=orig_fname))
write_cp2k(fhandle, "{}-supercell".format(pbase), supercell, tree)
for i, cell in enumerate(cells_with_displacements):
filename = "{pre_filename}-{0:0{width}}.inp".format(
suffix = "{pre_filename}-{0:0{width}}".format(
i + 1,
pre_filename=pre_filename,
width=width)
write_cp2k(filename, cell)
with open("{}-{}{}".format(fbase, suffix, fext), 'w') as fhandle:
fhandle.write("""\
# Generated by Phonopy, based on {fname}
# Merged configuration with displacements
""".format(fname=orig_fname))
write_cp2k(fhandle, "{}-{}".format(pbase, suffix), cell, tree)
def get_cp2k_structure(atoms):
"""Convert the atoms structure to a CP2K input file skeleton string"""
def write_cp2k(fhandle, project_name, atoms, tree):
"""Merge the new the atoms structure with the configuration tree to a new CP2K input file
from cp2k_tools.generator import dict2cp2k
:param fhandle: open file handle to which the routine will write to
:param project_name: the project name to use (CP2K uses that as prefix for generated files)i
:param atoms: the Atoms objects to use
:param tree: the configuration tree as returned from CP2KInputParser
"""
# CP2K's default unit is angstrom, convert it, but still declare it explictly:
cp2k_cell = {sym: ('[angstrom]',) + tuple(coords) for sym, coords in zip(('a', 'b', 'c'), atoms.get_cell()*Bohr)}
cp2k_cell['periodic'] = 'XYZ' # anything else does not make much sense
cp2k_coord = {
'scaled': True,
'*': [[sym] + list(coord) for sym, coord in zip(atoms.get_chemical_symbols(), atoms.get_scaled_positions())],
}
from cp2k_input_tools.generator import CP2KInputGenerator
generator = CP2KInputGenerator()
return dict2cp2k(
{
'global': {
'run_type': 'ENERGY_FORCE',
},
'force_eval': {
'subsys': {
'cell': cp2k_cell,
'coord': cp2k_coord,
},
'print': {
'forces': {
'filename': 'forces',
},
},
},
tree['+global']['run_type'] = 'ENERGY_FORCE'
tree['+global']['project_name'] = project_name
force_eval = tree['+force_eval'][0]
subsys = force_eval['+subsys']
# if the original input contained scaled positions, continue with scaled positions
if subsys['+coord'].get('scaled', False):
cp2k_coord = {
'scaled': True,
'*': ["{sym} {x} {y} {z}".format(sym=sym, x=coord[0], y=coord[1], z=coord[2])
for sym, coord in zip(atoms.get_chemical_symbols(), atoms.get_scaled_positions())],
}
)
# ... otherwise use absolute positions
else:
cp2k_coord = {
'*': ["{sym} {x} {y} {z}".format(sym=sym, x=coord[0], y=coord[1], z=coord[2])
for sym, coord in zip(atoms.get_chemical_symbols(), atoms.get_positions())],
}
subsys['+cell']['a'] = list(atoms.get_cell()[0])
subsys['+cell']['b'] = list(atoms.get_cell()[1])
subsys['+cell']['c'] = list(atoms.get_cell()[2])
subsys['+cell']['periodic'] = 'XYZ' # anything else does not make much sense
subsys['+coord'] = cp2k_coord # overwriting the coordinates
if not '+print' in force_eval:
force_eval['+print'] = {}
if not '+forces' in force_eval['+print']:
force_eval['+print']['+forces'] = {}
force_eval['+print']['+forces']['filename'] = "forces" # uses the project name as base with 'forces' as suffix
for line in generator.line_iter(tree):
fhandle.write("{line}\n".format(line=line))

View File

@ -66,7 +66,7 @@ AbinitToTHz = sqrt(EV/(AMU*Bohr))/Angstrom/(2*pi)/1e12 # [THz] 21.49068
PwscfToTHz = sqrt(Rydberg*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 108.97077
ElkToTHz = sqrt(Hartree*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 154.10794
SiestaToTHz = sqrt(EV/(AMU*Bohr))/Angstrom/(2*pi)/1e12 # [THz] 21.49068
CP2KToTHz = ElkToTHz # CP2K uses a.u. units (Hartree/Bohr)
CP2KToTHz = sqrt(Hartree*EV/(AMU*Bohr))/Angstrom/(2*pi)/1e12 # CP2K uses a.u. for forces but Angstrom for distances
CrystalToTHz = VaspToTHz
DftbpToTHz = sqrt(Hartree*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 154.10794344
dftbpToBohr = 0.188972598857892E+01

View File

@ -182,6 +182,7 @@ if __name__ == '__main__':
url='http://atztogo.github.io/phonopy/',
packages=packages_phonopy,
install_requires=['numpy', 'PyYAML', 'matplotlib', 'h5py'],
extras_require={'cp2k': ['cp2k-input-tools']},
provides=['phonopy'],
scripts=scripts_phonopy,
ext_modules=ext_modules_phonopy,