From 88c3b1212964a01f09156e8d858bba8ca1ca509d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tiziano=20M=C3=BCller?= Date: Mon, 16 Sep 2019 09:32:39 +0200 Subject: [PATCH] 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) --- example/Si-CP2K/README | 20 -- example/Si-CP2K/README.md | 26 +++ .../Si-CP2K/Si-supercell-001-forces-1_0.xyz | 130 +++++------ example/Si-CP2K/disp.yaml | 204 ---------------- phonopy/interface/calculator.py | 27 +-- phonopy/interface/cp2k.py | 221 ++++++++++++------ phonopy/units.py | 2 +- setup.py | 1 + 8 files changed, 251 insertions(+), 380 deletions(-) delete mode 100644 example/Si-CP2K/README create mode 100644 example/Si-CP2K/README.md delete mode 100644 example/Si-CP2K/disp.yaml diff --git a/example/Si-CP2K/README b/example/Si-CP2K/README deleted file mode 100644 index 6e87dfe0..00000000 --- a/example/Si-CP2K/README +++ /dev/null @@ -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" diff --git a/example/Si-CP2K/README.md b/example/Si-CP2K/README.md new file mode 100644 index 00000000..2286d3e9 --- /dev/null +++ b/example/Si-CP2K/README.md @@ -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" +``` diff --git a/example/Si-CP2K/Si-supercell-001-forces-1_0.xyz b/example/Si-CP2K/Si-supercell-001-forces-1_0.xyz index 1246be75..2f478a65 100644 --- a/example/Si-CP2K/Si-supercell-001-forces-1_0.xyz +++ b/example/Si-CP2K/Si-supercell-001-forces-1_0.xyz @@ -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 diff --git a/example/Si-CP2K/disp.yaml b/example/Si-CP2K/disp.yaml deleted file mode 100644 index ab87c460..00000000 --- a/example/Si-CP2K/disp.yaml +++ /dev/null @@ -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 \ No newline at end of file diff --git a/phonopy/interface/calculator.py b/phonopy/interface/calculator.py index 3b3c591c..01224cb7 100644 --- a/phonopy/interface/calculator.py +++ b/phonopy/interface/calculator.py @@ -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) diff --git a/phonopy/interface/cp2k.py b/phonopy/interface/cp2k.py index 862e98e2..412f349c 100644 --- a/phonopy/interface/cp2k.py +++ b/phonopy/interface/cp2k.py @@ -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)) diff --git a/phonopy/units.py b/phonopy/units.py index 10f54d2e..73c3cd3d 100644 --- a/phonopy/units.py +++ b/phonopy/units.py @@ -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 diff --git a/setup.py b/setup.py index f91c3a76..0c674bf3 100644 --- a/setup.py +++ b/setup.py @@ -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,