Solid-state Modeling Packages
To obtain the Hessian matrix from some solid-state modeling packages, the latest version of Phonopy package needs to be installed first.
Note
Although the following examples are on the calculation of solid systems with different solid-state modeling packages, this does not mean that molecular systems calculated by these packages cannot be analyzed by LModeA-nano.
In fact, molecular systems calculated by solid-state modeling packages can be analyzed in the same way as solids except that the Dimensions needs to be set to 0 in LModeA-nano.
VASP
Tested with VASP 5.4.4.
Method 1
The first step is cell optimization. For 3D solid system, an example INCAR file is
PREC = Accurate
GGA = MK
LUSE_VDW = .TRUE.
PARM1 = 0.1234
PARM2 = 0.711357
Zab_vdW = -1.8867
AGGAC = 0.000
LASPH = .TRUE.
IBRION = 1
NSW = 3000
NELMIN = 5
ISIF = 3
ENCUT = 1000.000000
EDIFF = 1.000000e-08
EDIFFG = -1.000000e-08
ISMEAR = 0
SIGMA = 0.1
IALGO = 38
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
NPAR = 4
After the cell optimization (ISIF=3) of a 3D solid structure, the vibrational analysis with VASP needs to be performed.
An example INCAR file for vibrational analysis is
PREC = Accurate
GGA = MK
LUSE_VDW = .TRUE.
PARM1 = 0.1234
PARM2 = 0.711357
Zab_vdW = -1.8867
AGGAC = 0.000
LASPH = .TRUE.
IBRION = 5
NFREE=2
POTIM=0.005
NELMIN = 5
ISIF = 3
ENCUT = 1000.000000
EDIFF = 1.000000e-08
EDIFFG = -1.000000e-08
ISMEAR = 0
SIGMA = 1.000000e-02
IALGO = 38
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
The most critical parameters are
IBRION=5- finite difference method to calculate the Hessian matrix without considering symmetryNFREE=2- use central differencesPOTIM=0.005- step size for finite difference in Å
When the vibrational analysis calculation is done, check the lowest frequencies by
$ grep cm OUTCAR | tail -n 6
115 f = 0.899811 THz 5.653679 2PiTHz 30.014464 cm-1 3.721321 meV
116 f = 0.729527 THz 4.583752 2PiTHz 24.334394 cm-1 3.017081 meV
117 f = 0.483556 THz 3.038272 2PiTHz 16.129693 cm-1 1.999828 meV
118 f/i= 0.013493 THz 0.084780 2PiTHz 0.450086 cm-1 0.055804 meV
119 f/i= 0.036011 THz 0.226262 2PiTHz 1.201188 cm-1 0.148928 meV
120 f/i= 0.056606 THz 0.355667 2PiTHz 1.888177 cm-1 0.234104 meV
Such result shows that the optimization is successful because the last three frequencies (close to zero cm \(^{-1}\)) correspond to the translation modes of the whole primitive unit cell and the remaining frequencies are all positive.
In the same directory of OUTCAR, another output file vasprun.xml is also generated. Run:
$ phonopy --fc vasprun.xml
...
FORCE_CONSTANTS has been created from vasprun.xml.
...
Two files are needed for later use
POSCAR- input structure for vibrational analysis calculationFORCE_CONSTANTS- Hessian matrix extracted by Phonopy from .xml file
Compose an input file vasp.inp for LModeA-nano
@VASP
POSCON = POSCAR
FCPHONOPY = FORCE_CONSTANTS
Following these steps
open a new PyMOL window and launch LModeA-nano by clicking Plugin → LModeA-nano
click the … (Browse) button and select the
vasp.inpfileclick Load
change Dimensions to
3and click Confirm buttonchange the supercell size percentage from
110%to140%or other values, and click Partial Supercell button
Then local mode analysis can be performed as described in quickstart examples.
Method 2 - For large cells
For primitive unit cells having more than 70 atoms, the calculation of Hessian matrix by finite difference method can be time-consuming because it requires 6N force evaluation steps (N is number of atoms in the cell). The above Method 1 performs all force evaluations on a single computational node. If there is a computational cluster, the following method is faster because it may distribute the force evaluations to different nodes.
After the cell optimization (ISIF=3) of a 3D solid structure, save the CONTCAR file as POSCAR in a new directory.
Run Phonopy
$ phonopy -d --dim="1 1 1" --pm --amplitude 0.005
...
Creating displacements
Plus Minus displacement: full plus minus directions
Settings:
Supercell: [1 1 1]
...
A series of displaced unit cell structure are generated: POSCAR-001, POSCAR-002, .... Each structrue will have its forces evaluated in a single VASP calculation.
The INCAR file for force evaluation can look like
PREC = Accurate
GGA = MK
LUSE_VDW = .TRUE.
PARM1 = 0.1234
PARM2 = 0.711357
Zab_vdW = -1.8867
AGGAC = 0.000
LASPH = .TRUE.
IBRION = -1
NELMIN = 5
ISIF = 3
ENCUT = 1000.000000
EDIFF = 1.000000e-08
EDIFFG = -1.000000e-08
ISMEAR = 0
SIGMA = 1.000000e-02
IALGO = 38
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
It is recommended to create a separate directory for each structure
for id in POSCAR-*
do
mkdir disp-${id}
cp ${id} disp-${id}/POSCAR
cp POTCAR KPOINTS INCAR disp-${id}/
done
When submitting the VASP calculations, each computational node can handle a few structures, like
for i in $(seq -f "%03g" 1 100)
do
echo disp-POSCAR-${i}
cd disp-POSCAR-${i}
vasp_std > vasp-${i}.out
cd ..
done
seq -f "%03g" 1 100 means to loop over 001 to 100.
After the force evaluation of all displaced structures are finished, create FORCE_SET using Phonopy
$ phonopy -f disp-POSCAR-{001..xxx}/vasprun.xml
xxx is the total number of displaced structures
Then use Phonopy to create FORCE_CONSTANTS file
$ phonopy --dim="1 1 1" --writefc
With POSCAR and FORCE_CONSTANTS files, the local mode analysis can be performed with LModeA-nano like in Method 1.
CP2K
Tested with CP2K 8.2. The input file template can be generated by Multiwfn package.
Method 1
The input file for cell optimization of 3D solid is like
#Generated by Multiwfn
&GLOBAL
PROJECT k2reh9-refine
PRINT_LEVEL LOW
RUN_TYPE CELL_OPT
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&SUBSYS
&CELL
A 9.691 0.000 0.000
B -4.810 8.413 0.000
C 0.000 0.000 5.706
PERIODIC XYZ #Direction of applied PBC (geometry aspect)
&END CELL
&COORD
K -1.2558488724 2.2006561653 2.8528315416
K 3.6056710442 6.2147826995 2.8528315383
K 2.5337691322 0.0019792003 2.8528315571
K -2.8321552830 4.9560287939 0.0000000311
K 2.0066157818 3.4586295569 -0.0000000269
K 5.7081521839 0.0010789324 0.0000000278
Re 4.8571785554 2.8039248183 2.8528315499
Re 0.0234382187 5.6083829765 2.8528315526
Re -0.0000638786 -0.0001099606 -0.0000000055
H 6.0270652495 2.5945198607 4.0766192575
H 4.4539534701 3.9229323799 4.0757259845
H 4.0894727226 1.8947602744 4.0751317071
H -0.3848103609 4.4906671571 1.6305327076
H 4.0894727162 1.8947602632 1.6305314009
H 6.0270652713 2.5945198686 1.6290438503
H 4.4539534568 3.9229323662 1.6299370985
H -0.7389977868 6.5200724154 4.0766186544
H 1.1950192311 5.8137176574 4.0757278734
H 1.6781788362 -0.0039857278 -0.0000000099
H 4.0388741329 6.9614774107 -0.0000000006
H -0.8363820802 1.4549103043 -0.0000000041
H 8.4818032881 0.0022916350 1.1978086117
H 0.6071583592 1.0465002108 4.5092456422
H -4.2077347386 7.3645633109 4.5078562001
H -4.2077347299 7.3645632719 1.1978068819
H 0.6071583649 1.0465001732 1.1964174406
H 8.4818032732 0.0022916783 4.5078544943
H -1.6284528162 5.3254766136 2.8528315589
H 0.6034005602 7.1805921593 2.8528315404
H 1.0952417964 4.3204237909 2.8528315625
H 3.2070901557 3.0951493757 2.8528315485
H 5.9342251097 4.0877344629 2.8528315420
H 5.4314516435 1.2294438796 2.8528315643
H -0.3848103509 4.4906671776 4.0751304240
H -0.7389977923 6.5200723974 1.6290444364
H 1.1950192243 5.8137176406 1.6299352197
&END COORD
&KIND K
ELEMENT K
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND Re
ELEMENT Re
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
# WFN_RESTART_FILE_NAME k2reh9-opt-RESTART.wfn
CHARGE 0 #Net charge
MULTIPLICITY 1 #Spin multiplicity
&KPOINTS
SCHEME MONKHORST-PACK 5 5 8
SYMMETRY F #If using symmetry to reduce the number of k-points
&END KPOINTS
&QS
EPS_DEFAULT 1E-10 #This is default. Set all EPS_xxx to values such that the energy will be correct up to this value
&END QS
&POISSON
PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
PSOLVER PERIODIC #The way to solve Poisson equation
&END POISSON
&XC
&XC_FUNCTIONAL PBE
&PBE
PARAMETRIZATION REVPBE
&END PBE
&END XC_FUNCTIONAL
&END XC
&MGRID
CUTOFF 825
REL_CUTOFF 75
&END MGRID
&SCF
MAX_SCF 128
EPS_SCF 1.0E-06 #Convergence threshold of density matrix during SCF
# SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
&DIAGONALIZATION
ALGORITHM STANDARD #Algorithm for diagonalization. DAVIDSON is faster for large systems
&END DIAGONALIZATION
&MIXING #How to mix old and new density matrices
METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
&END MIXING
&PRINT
&RESTART #Use "&RESTART OFF" can prevent generating wfn file
BACKUP_COPIES 0 #Maximum number of backup copies of wfn file
&END RESTART
&END PRINT
&END SCF
&END DFT
STRESS_TENSOR ANALYTICAL
&END FORCE_EVAL
&MOTION
&CELL_OPT
MAX_ITER 250 #Maximum number of geometry optimization
EXTERNAL_PRESSURE 1.01325 #External pressure for cell optimization (bar)
CONSTRAINT NONE #Can be e.g. Z, XY to fix corresponding cell length
KEEP_ANGLES F #If T, then cell angles will be kepted
KEEP_SYMMETRY F #If T, crystal symmetry will be kepted, and symmetry should be specified in &CELL
TYPE DIRECT_CELL_OPT #Geometry and cell are optimized at the same time. Can also be GEO_OPT, MD
#The following thresholds of optimization convergence are the default ones
MAX_DR 0.000060 #Maximum geometry change
RMS_DR 0.000040 #RMS geometry change
MAX_FORCE 0.000015 #Maximum force
RMS_FORCE 0.000010 #RMS force
PRESSURE_TOLERANCE 100 #Pressure tolerance (w.r.t EXTERNAL_PRESSURE)
OPTIMIZER BFGS #Can also be CG (more robust for difficult cases) or LBFGS
&BFGS
TRUST_RADIUS 0.2 #Trust radius (maximum stepsize) in Angstrom
# RESTART_HESSIAN T #If read initial Hessian, uncomment this line and specify the file in the next line
# RESTART_FILE_NAME to_be_specified
&END BFGS
&END CELL_OPT
&PRINT
&RESTART
BACKUP_COPIES 0 #Maximum number of backing up restart file
&END RESTART
&RESTART_HISTORY
&EACH
CELL_OPT 0 #How often a history .restart file is generated, 0 means never
&END EACH
&END RESTART_HISTORY
&CELL
&END CELL
&END PRINT
&END MOTION
With the optimized cell structure, the vibrational analysis is required
#Generated by Multiwfn
&GLOBAL
PROJECT k2-refine-freq
PRINT_LEVEL MEDIUM
RUN_TYPE VIBRATIONAL_ANALYSIS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&SUBSYS
&CELL
A 9.7154601050 0.0000000000 0.0000000000
B -4.8583032251 8.4134784197 0.0000000000
C -0.0000048135 -0.0000078762 5.6945465926
PERIODIC XYZ #Direction of applied PBC (geometry aspect)
&END CELL
&COORD
K -1.2717444316 2.2034502965 2.8470104286
K 3.5860479359 6.2116264326 2.8470414297
K 2.5441673793 0.0005069845 2.8470105500
K -2.8664982103 4.9640103204 -0.0002691372
K 1.9919809386 3.4504437696 -0.0002481916
K 5.7322486247 -0.0000606935 -0.0002691208
Re 4.8578753089 2.8054541385 2.8470302019
Re 0.0003345425 5.6097335267 2.8470301311
Re 0.0003406364 0.0005363186 -0.0002407693
H 6.0318408651 2.6225498479 4.0708725666
H 4.4291201342 3.9136919383 4.0706909878
H 4.1125879533 1.8801461820 4.0708571241
H -0.4282825066 4.5016125249 1.6232100820
H 4.1125863706 1.8801513315 1.6232012617
H 6.0318388962 2.6225540649 1.6231796732
H 4.4291142103 3.9136997301 1.6233758117
H -0.7451026988 6.5349089400 4.0708601931
H 1.1744732842 5.7926227251 4.0707068918
H 1.6786597913 0.0005874087 -0.0002624812
H 4.0184027906 6.9605294583 -0.0002696641
H -0.8388731978 1.4539493631 -0.0002624278
H 8.5073111255 0.0005134857 1.1978276666
H 0.6045199441 1.0470807705 4.4961739186
H -4.2536730257 7.3674658865 4.4962388277
H -4.2536811640 7.3674669181 1.1978123072
H 0.6045200337 1.0470806825 1.1978988291
H 8.5073057336 0.0005201478 4.4962231782
H -1.6564785794 5.3554871447 2.8470329264
H 0.6084571955 7.1717438292 2.8470366227
H 1.0487876716 4.3019244622 2.8470389912
H 3.2010232060 3.0594336945 2.8470390703
H 5.9064519076 4.1131982559 2.8470367529
H 5.4661799457 1.2435429741 2.8470328975
H -0.4282878791 4.5016116953 4.0708480056
H -0.7450971903 6.5349082052 1.6231930653
H 1.1744841476 5.7926214329 1.6233587569
&END COORD
&KIND K
ELEMENT K
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND Re
ELEMENT Re
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
# WFN_RESTART_FILE_NAME k2-freq-template-RESTART.wfn
CHARGE 0 #Net charge
MULTIPLICITY 1 #Spin multiplicity
&KPOINTS
SCHEME MONKHORST-PACK 5 5 8
SYMMETRY F #If using symmetry to reduce the number of k-points
&END KPOINTS
&QS
EPS_DEFAULT 1E-10 #This is default. Set all EPS_xxx to values such that the energy will be correct up to this value
&END QS
&POISSON
PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
PSOLVER PERIODIC #The way to solve Poisson equation
&END POISSON
&XC
&XC_FUNCTIONAL PBE
&PBE
PARAMETRIZATION REVPBE
&END PBE
&END XC_FUNCTIONAL
&END XC
&MGRID
CUTOFF 825
REL_CUTOFF 75
&END MGRID
&SCF
MAX_SCF 128
EPS_SCF 1.0E-07 #Convergence threshold of density matrix during SCF
# SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
&DIAGONALIZATION
ALGORITHM STANDARD #Algorithm for diagonalization. DAVIDSON is faster for large systems
&END DIAGONALIZATION
&MIXING #How to mix old and new density matrices
METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
&END MIXING
&END SCF
&PRINT
&MOMENTS
&END MOMENTS
&END PRINT
&END DFT
&END FORCE_EVAL
&MOTION
&PRINT
&RESTART
BACKUP_COPIES 0 #Maximum number of backing up restart file
&END RESTART
&END PRINT
&END MOTION
&VIBRATIONAL_ANALYSIS
DX 0.008 #Step size of finite difference. This is default (Bohr)
NPROC_REP 1 #Number of processors to be used per replica. This is default
TC_PRESSURE 101325 #1 atm. Pressure for calculate thermodynamic data (Pa)
TC_TEMPERATURE 298.15 #Temperature for calculate thermodynamic data (K)
THERMOCHEMISTRY #Print thermochemistry information (only correct for gas molecule!)
INTENSITIES F #Calculate IR intensities
&PRINT
&HESSIAN high
&END HESSIAN
&MOLDEN_VIB #Output .mol (Molden file) for visualization vibrational modes
&END MOLDEN_VIB
&END PRINT
&END VIBRATIONAL_ANALYSIS
Critical parameters include
DX 0.008- this specifies the step size of finite differenceHESSIAN high- this prints the Hessian matrix in the CP2K output file
After the vibrational analysis is done, check the vibrational frequency values to make sure the structure is at a local minimum energetically.
Compose an input file cp2k.inp for LModeA-nano
@cp2k
mode = native
OUTPUT = k2reh9-refine-freq.log
INPUT = k2reh9-refine-freq.inp
The INPUT and OUTPUT specify the input and output files of vibrational analysis calculation.
Follow these steps to perform local mode analysis with LModeA-nano
open a new PyMOL window and launch LModeA-nano by clicking Plugin → LModeA-nano
click the … (Browse) button and select the
cp2k.inpfilechange the program drop-down list from VASP to CP2K
click Load
change Dimensions to
3and click Confirm buttonchange the supercell size percentage from
110%to140%or other values, and click Partial Supercell button
Then local mode analysis can be performed as described in quickstart examples.
Method 2 - For large cells
Similar to Method 2 with VASP, the general idea of this method is ask Phonopy package to generate the CP2K input files for the displaced structures starting from an already optimized structure. Then force evaluations are performed on these displaced structures in order to obtain Hessian.
A template CP2K input file having the optimized cell structure is needed, and the format of this input should look like the following. Otherwise, Phonopy package cannot parse it correctly. More examples of such CP2K input can be found in the Phonopy repo on GitHub: NaCl and Si.
! whether to do kpoints or not
! without kpoints it will be significantly faster,
! but some crossings in the band structure will be missing
@SET WITH_KP yes
! set the following to yes to do a cell optimization instead of energy/force calculation
! do NOT enable this when using this input file as template for phonopy!
@SET DO_CELLOPT no
&GLOBAL
@IF $DO_CELLOPT == yes
PROJECT solid-cellopt
RUN_TYPE CELL_OPT
@ENDIF
@IF $DO_CELLOPT == no
PROJECT na2reh9-phonopy
RUN_TYPE ENERGY_FORCE
@ENDIF
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
&QS
EPS_DEFAULT 1e-10
! We are using GAPW instead of the default GPW to prevent numerical
! noise due to the finite grid size in the integration schemes in CP2K
! together with the small displacements used to determine the forces.
! Alternatively one could increase the CUTOFF to >5000 (depends on basis set)
! to ensure that the higher exponents in the basis sets are well represented on the grid.
&END QS
&POISSON
PERIODIC XYZ
PSOLVER PERIODIC
&END POISSON
&SCF
EPS_SCF 1e-07
MAX_SCF 80
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING #How to mix old and new density matrices
METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
&END MIXING
&END SCF
@IF $WITH_KP == yes
&KPOINTS
SCHEME MONKHORST-PACK 5 5 8
SYMMETRY .FALSE.
&END KPOINTS
@ENDIF
&XC
&XC_FUNCTIONAL PBE
&PBE
PARAMETRIZATION REVPBE
&END PBE
&END XC_FUNCTIONAL
&END XC
&MGRID
REL_CUTOFF 80
CUTOFF 900
&END MGRID
&END DFT
&SUBSYS
&CELL
A 8.8901326513 0.0000000000 0.0000000000
B -4.4739104629 7.6824375325 0.0000000000
C 0.0000001000 0.0000003943 5.2483716421
PERIODIC XYZ
&END CELL
&COORD
Na -1.1859069865 2.0379623261 2.6241869010
Na 3.2449119045 5.6448521398 2.6241867646
Na 2.3578700048 0.0007892838 2.6241868726
Na -2.6256486318 4.5068346666 0.0000001680
Na 1.8251693178 3.1750730554 0.0000001576
Na 5.2158593380 -0.0009231392 0.0000001673
Re 4.4364564557 2.5642868258 2.6242178317
Re -0.0167453884 5.1242663051 2.6242177424
Re 0.0017170475 0.0030111053 0.0000094252
H 5.6222385887 2.4929720119 3.8444638123
H 3.9052360919 3.6232613245 3.8472931155
H 3.7846319842 1.5754147357 3.8473354516
H -0.5432325121 4.0633761049 1.4010781033
H 3.7846337784 1.5754183851 1.4010976360
H 5.6222353787 2.4929729080 1.4039700228
H 3.9052378288 3.6232592016 1.4011411061
H -0.6751167153 6.1130464337 3.8444723455
H 1.1657297010 5.1981428489 3.8472647704
H 1.6763457738 0.0076515368 0.0000010890
H 3.5835080375 6.2338251923 0.0000016665
H -0.8370053879 1.4524719986 0.0000010868
H 7.6815358956 0.0008872303 1.1958253823
H 0.6038026414 1.0503904668 4.0506195218
H -3.8649471720 6.6385181149 4.0525569576
H -3.8649468269 6.6385180090 1.1958373681
H 0.6038029996 1.0503908521 1.1977742281
H 7.6815360048 0.0008870488 4.0525690768
H -1.6822864070 5.0230528750 2.6242153250
H 0.7327983274 6.6152712210 2.6242163261
H 0.9014162231 3.7301543776 2.6242158401
H 2.7696859375 2.6561823231 2.6242158698
H 5.3476415856 3.9623777672 2.6242164032
H 5.1871788616 1.0740919399 2.6242154475
H -0.5432348504 4.0633727533 3.8473548145
H -0.6751143293 6.1130441536 1.4039612322
H 1.1657270304 5.1981432445 1.4011693676
&END COORD
&KIND Na
ELEMENT Na
BASIS_SET DZVP-MOLOPT-SR-GTH # use an AE basis optimized for solids, any other (matching the pseudo) will work too, though
POTENTIAL GTH-PBE # one could also use a pseudopotentials (non-AE)
! while not strictly required, this should make the initial guess more accurate (Na+ Cl-)
&END KIND
&KIND Re
ELEMENT Re
BASIS_SET DZVP-MOLOPT-SR-GTH # use an AE basis optimized for solids, any other (matching the pseudo) will work too, though
POTENTIAL GTH-PBE # one could also use a pseudopotentials (non-AE)
! while not strictly required, this should make the initial guess more accurate (Na+ Cl-)
&END KIND
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH # use an AE basis optimized for solids, any other (matching the pseudo) will work too, though
POTENTIAL GTH-PBE # one could also use a pseudopotentials (non-AE)
! while not strictly required, this should make the initial guess more accurate (Na+ Cl-)
&END KIND
&END SUBSYS
STRESS_TENSOR ANALYTICAL
@IF $DO_CELLOPT == no
&PRINT
&STRESS_TENSOR
ADD_LAST NUMERIC
FILENAME stress_tensor
&END STRESS_TENSOR
&FORCES
ADD_LAST NUMERIC
FILENAME forces
&END FORCES
&END PRINT
@ENDIF
&END FORCE_EVAL
@IF $DO_CELLOPT == yes
&MOTION
&CELL_OPT
KEEP_ANGLES .TRUE.
MAX_FORCE 1.0E-10
&END CELL_OPT
&END MOTION
@ENDIF
Then run
$ phonopy --cp2k -d --dim="1 1 1" --pm --amplitude 0.005 -c na2reh9-template.inp
The parameter after -c is the above-mentioned CP2K template input file containing the optimized structure. Phonopy generates a series of CP2K input files for force evaluation: na2reh9-template-supercell-001.inp, na2reh9-template-supercell-002.inp, .... This step will also generate a phonopy_disp.yaml file.
The generated CP2K input can look like
# Generated by Phonopy, based on na2reh9-template.inp
# Merged configuration with displacements
&FORCE_EVAL
METHOD QUICKSTEP
STRESS_TENSOR ANALYTICAL
&PRINT
&FORCES
ADD_LAST NUMERIC
FILENAME forces
&END FORCES
&STRESS_TENSOR
ADD_LAST NUMERIC
FILENAME stress_tensor
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND Re
ELEMENT Re
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND Na
ELEMENT Na
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&COORD
Na -1.1809069865 2.0379623261 2.624186901
Na 3.2449119044999986 5.6448521398 2.6241867646
Na 2.3578700047999996 0.0007892838 2.6241868725999993
Na 6.2644840195 4.5068346666 1.6799999999999997e-07
Na 1.825169317800001 3.1750730554 1.576e-07
Na 0.7419488750999992 7.6815143933 1.673e-07
Re 4.436456455700002 2.5642868258 2.6242178317
Re -0.01674538839999951 5.124266305100002 2.6242177424
Re 0.0017170474999999995 0.003011105300000001 9.425199999999999e-06
H 5.6222385887 2.4929720119 3.844463812299999
H 3.9052360918999995 3.6232613245 3.847293115499999
H 3.7846319841999994 1.5754147357000003 3.847335451599999
H -0.5432325121 4.0633761049 1.4010781033
H 3.7846337784 1.5754183851 1.401097636
H 5.6222353787 2.492972907999999 1.4039700228
H 3.9052378287999994 3.6232592016 1.4011411060999996
H -0.6751167152999991 6.113046433700002 3.8444723455
H 1.165729701000001 5.198142848899998 3.8472647704
H 1.6763457738 0.0076515368000000025 1.089e-06
H 3.583508037500002 6.2338251923 1.6665e-06
H -0.8370053878999996 1.4524719986 1.0867999999999997e-06
H 7.6815358956 0.0008872303000000002 1.1958253823
H 0.6038026414 1.0503904668 4.0506195218
H -3.864947171999999 6.6385181149 4.0525569576
H -3.864946826899999 6.638518009 1.1958373681
H 0.6038029995999996 1.0503908521 1.1977742281
H 7.6815360048 0.0008870488 4.052569076799999
H -1.682286407 5.023052875 2.6242153249999998
H 0.7327983273999987 6.615271221 2.624216326099999
H 0.9014162231000011 3.730154377600001 2.6242158401
H 2.7696859375000003 2.6561823231 2.6242158698
H 5.3476415856 3.9623777672 2.6242164032
H 5.1871788616 1.0740919399 2.6242154474999997
H -0.5432348503999994 4.0633727533 3.8473548145
H -0.6751143292999994 6.1130441536 1.4039612321999997
H 1.1657270304000014 5.1981432445 1.4011693676
&END COORD
&CELL
A 8.8901326513 0.0 0.0
B -4.4739104629 7.6824375325 0.0
C 1e-07 3.943e-07 5.2483716421
PERIODIC XYZ
&END CELL
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
REL_CUTOFF 80.0
CUTOFF 900.0
&END MGRID
&XC
&XC_FUNCTIONAL PBE
&PBE
PARAMETRIZATION REVPBE
&END PBE
&END XC_FUNCTIONAL
&END XC
&KPOINTS
SCHEME MONKHORST-PACK 5 5 8
SYMMETRY .FALSE.
&END KPOINTS
&SCF
EPS_SCF 1e-07
MAX_SCF 80
&MIXING
METHOD BROYDEN_MIXING
ALPHA 0.4
NBUFFER 8
&END MIXING
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&END SCF
&POISSON
PERIODIC XYZ
POISSON_SOLVER PERIODIC
&END POISSON
&QS
EPS_DEFAULT 1e-10
&END QS
&END DFT
&END FORCE_EVAL
&GLOBAL
PROJECT_NAME na2reh9-phonopy-supercell-001
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL
After the force evaluation for these generated input files is completed, it creates a series files containing the force information: na2reh9-phonopy-supercell-001-forces-1_0.xyz, na2reh9-phonopy-supercell-002-forces-1_0.xyz, ...
Make sure that these .xyz files are collected in one directory and this directory should contain phonopy_disp.yaml file.
Run Phonopy to create the FORCE_SETS file
$ phonopy --cp2k -f na2reh9-phonopy-supercell-{001..xxx}-forces-1_0.xyz
xxx is the total number of displaced structures generated by Phonopy
Then use Phonopy to generate the FORCE_CONSTANTS file
$ phonopy --dim="1 1 1" --writefc
Compose the input file cp2k.inp for LModeA-nano
@cp2k
mode = phonopy
FCPHONOPY = FORCE_CONSTANTS
INPUT = na2reh9-template.inp
The INPUT parameter specifies the template CP2K input file containing the optimized structure.
Quantum ESPRESSO
Tested with QE 6.4.1.
There are two ways to do local mode analysis on the results from QE. The first one is to utilize the density functional perturbation theory (DFPT) implemented in QE to obtain the force contant matrix at Gamma point (q=0). The other method known as “frozon phonon” is to use finite difference approach to generate a series of force sets in order to construct the Hessian matrix by using Phonopy package.
These two approaches are equivalent in principle. While the DFPT method in QE has limited support for levels of theory (e.g. hybrid density functional), one has to use the finite difference approach if DFPT is not applicable. Therefore, it is recommended to employ the frozon phonon method.
Suppose that users have obtained the optimized primitive cell structure of a solid, the following tutorial skips this cell optimization step.
DFPT
This method includes three steps
SCF calculation with
pw.xDFPT calculation with
ph.xgenerate force constant matrix with
q2r.x
Based on an optimized primitive unit cell structure, prepare an SCF input file diamond-scf.inp
&CONTROL
pseudo_dir = '.'
outdir = './outp'
calculation = 'scf'
prefix='carbon'
/
&SYSTEM
ibrav = 0
nat = 2
ntyp = 1
ecutwfc = 40.0
ecutrho = 400.0
/
&ELECTRONS
scf_must_converge = .true.
/
ATOMIC_SPECIES
C 12.0107 C.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (angstrom)
2.522557697 0.000000000 -0.000000000
1.261278849 2.184599048 -0.000000000
1.261278849 0.728199683 2.060840365
ATOMIC_POSITIONS (angstrom)
C 3.783836546 2.184599048 1.545801941
C 2.522557698 1.456399366 1.030248516
K_POINTS automatic
5 5 5 1 1 1
The CELL_PARAMETERS and ATOMIC_POSITIONS have to be in the unit of angstrom.
Run DFPT calculations. The program to run this job is PHonon (ph.x). The input file needs to specify the same outdir and prefix as the SCF job did. In addition, specify the Gamma point (q=[0,0,0]) in the end of the input file.
phonons of diamond at Gamma
&inputph
tr2_ph=1.0d-14,
epsil=.true.,
prefix='carbon',
fildyn='dynmat.dyn1',
amass(1)=12.0107,
outdir='./outp'
/
0.0 0.0 0.0
The fildyn specifies the name of the dynamical matrix file as an output.
Create a text file named dynmat.dyn0
1 1 1
1
Run
$ cp dynmat.dyn1 dynmat.dyn
Prepare the input file for q2r.x
&input
fildyn='dynmat.dyn',
zasr='simple',
flfrc='diamond.fc'
/
In this input file, the zasr controls the the acoustic sum rules. For solid systems, we need to set zasr to 'simple' or 'crystal'. For one-dimensional polymer systems, we need to set zasr to 'one-dim'. For isolated molecular systems, we need to set zasr to 'zero-dim'.
The parameter flfrc specifies the name of the interatomic force constant file generated from the calculation of q2r.x.
So far, we have obtained necessary data files for local mode analysis
diamond-scf.inp- input file for SCF calculationdiamond.fc- the force constant file specified byflfrc
Compose an input file qe.inp for local mode analysis
@QE
mode = dfpt
GeomIn = diamond-scf.in
q2rFC = diamond.fc
Follow these steps to perform local mode analysis with LModeA-nano
open a new PyMOL window and launch LModeA-nano by clicking Plugin → LModeA-nano
click the … (Browse) button and select the
qe.inpfilechange the program drop-down list from VASP to Quantum ESPRESSO
click Load
change Dimensions to
3and click Confirm buttonchange the supercell size percentage from
110%to140%or other values, and click Partial Supercell button
Then local mode analysis can be performed as described in quickstart examples.
Frozon Phonon
Given an already optimized cell structure, prepare a template QE input file diamond-temp.in
&CONTROL
pseudo_dir = '.'
calculation = 'scf'
tprnfor = .true.
tstress = .true.
/
&SYSTEM
ibrav = 0
nat = 2
ntyp = 1
ecutwfc = 40.0
ecutrho = 400.0
/
&ELECTRONS
scf_must_converge = .true.
/
ATOMIC_SPECIES
C 12.0107 C.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (angstrom)
2.522557697 0.000000000 -0.000000000
1.261278849 2.184599048 -0.000000000
1.261278849 0.728199683 2.060840365
ATOMIC_POSITIONS (crystal)
C 0.749972224 0.749972284 0.750083327
C 0.500027776 0.500027716 0.499916703
K_POINTS automatic
5 5 5 1 1 1
It is important to set tprnfor and tstress as force evaluation is desired to construct Hessian. In addition, the atomic positions have to be in the crystal format using the fractional coordinates as required by Phonopy later.
Then save upper and lower part as two separate files
diamond-temp-upper.in :
&CONTROL
pseudo_dir = '.'
calculation = 'scf'
tprnfor = .true.
tstress = .true.
/
&SYSTEM
ibrav = 0
nat = 2
ntyp = 1
ecutwfc = 40.0
ecutrho = 400.0
/
&ELECTRONS
scf_must_converge = .true.
/
diamond-temp-lower.in :
K_POINTS automatic
5 5 5 1 1 1
Use Phonopy to generate displaced structures
$ phonopy --qe -c diamond-temp.in -d --dim="1 1 1"
This will generates a series of files containing the displaced structures supercell-001.in, supercell-002.in, .... But these files only contains the structure information. Create complete QE input files with the following script
prefix="supercell-"
for a in supercell-*.in
do
b=${a#$prefix}
cat diamond-temp-upper.in $a diamond-temp-lower.in > diamond-${b}
done
This script can generate QE input files diamond-001.in, diamond-002.in, .... These input files differ from diamond-temp.in only in the structure information.
Run QE pw.x calculations with these generated input files, the output files are expected to be diamond-001.out, diamond-002.out, ...
Use Phonopy to generate the FORCE_SETS file
$ phonopy --qe -f diamond-{001..xxx}.out
xxx is the total number of displaced structures.
Generate FORCE_CONSTANTS file with Phonopy
$ phonopy --writefc
Compose a LModeA-nano input file qe.inp
@QE
mode = phonopy
GeomIn = diamond-temp.in
FCPHONOPY = FORCE_CONSTANTS
The remaining step to perform local mode analysis is the same as in DFPT method.
CASTEP
Tested with CASTEP 20.11. CASTEP is able to perform phonon (vibrational analysis) calculations with either DFPT or finite difference method.
Given the already optimized solid structure, prepare the .cell input for CASTEP
%BLOCK lattice_cart
Ang
2.5333173 0.0000000 0.0000000
1.2666587 2.1939171 0.0000000
1.2666587 0.7313057 2.0684936
%ENDBLOCK lattice_cart
%BLOCK positions_frac
C 0.749996 0.749996 0.750013
C 0.500004 0.500004 0.499987
%ENDBLOCK positions_frac
!
! Choose which pseudopotentials to use
! Either specify external files, or omit to generate a pseudopotential
!
%BLOCK species_pot
C NCP
%ENDBLOCK species_pot
!
! Analyse structure to determine symmetry
!
!symmetry_generate
!
! Specify M-P grid dimensions for electron wavevectors (K-points)
!
kpoint_mp_grid 5 5 5
!
! Offset the grid in the xy plane to include (0,0,0)
%block PHONON_KPOINT_LIST
0.0 0.0 0.0 1.0 ! Wavevector of phonon(s) to compute ( qx qy qz, weight)
%endblock PHONON_KPOINT_LIST
Then prepare the .param file for either DFPT
IPRINT 3
!
! Task keyword specifies type of calculation
task : phonon
!
! Treat the system as an insulator
!
fix_occupancy : true
!
! Tune for speed rather than memory conservation.
!
opt_strategy : speed
!
! Don't write unnecessary wavefunction files
!
num_dump_cycles : 0
!
! Choose exchange-correlation functional
!
xc_functional : PBE
!
! Give plane-wave cutoff. Also allowed keywords (COARSE,MEDIUM,FINE,PRECISE)
!
BASIS_PRECISION : FINE
!
! Choose size of FFT grid. 2.0 is fully converged. 1.75 is good enough for LDA
!
grid_scale 2.0
!
! Choose which SCF solver to use ("none" for allbands, "dm" for Density Mixing)
!
elec_method : dm
mixing_scheme : Pulay
mix_charge_amp : 0.6
!
! Convergence tolerance criterion - energy per atom.
!
elec_energy_tol : 1.0e-9 eV
!
! Turn on calculation of stress (off by default)
!
calculate_stress : true
phonon_sum_rule : true
or finite difference calculation
IPRINT 2
!
! Task keyword specifies type of calculation
task : phonon
PHONON_METHOD : FINITEDISPLACEMENT
!
! Treat the system as an insulator
!
fix_occupancy : true
!
! Tune for speed rather than memory conservation.
!
opt_strategy : speed
!
! Don't write unnecessary wavefunction files
!
num_dump_cycles : 0
!
! Choose exchange-correlation functional
!
xc_functional : PBE
!
! Give plane-wave cutoff. Also allowed keywords (COARSE,MEDIUM,FINE,PRECISE)
!
BASIS_PRECISION : FINE
!
! Choose size of FFT grid. 2.0 is fully converged. 1.75 is good enough for LDA
!
grid_scale 2.0
!
! Choose which SCF solver to use ("none" for allbands, "dm" for Density Mixing)
!
elec_method : dm
mixing_scheme : Pulay
mix_charge_amp : 0.6
!
! Convergence tolerance criterion - energy per atom.
!
elec_energy_tol : 1.0e-9 eV
!
! Turn on calculation of stress (off by default)
!
calculate_stress : true
phonon_sum_rule : true
The major output is .castep file. Then compose an input file castep.inp for LModeA-nano.
@castep
output = xxxx.castep
Follow these steps
open a new PyMOL window and launch LModeA-nano by clicking Plugin → LModeA-nano
click the … (Browse) button and select the
qe.inpfilechange the program drop-down list from VASP to CASTEP
click Load
change Dimensions to
3and click Confirm buttonchange the supercell size percentage from
110%to140%or other values, and click Partial Supercell button
Then local mode analysis can be performed as described in quickstart examples.
CRYSTAL
Tested with CRYSTAL17.
To optimize the primitive cell structure, the input file can be
graphene
SLAB
1
2.456 2.456 120.00
2
6 0.000000000 0.000000000 0.000000000
6 0.333333333 0.666666666 0.000000000
OPTGEOM
TOLDEG
0.0000050
TOLDEX
0.000010
TOLDEE
1E-9
END
END
6 6
0 0 5 2 1
1238.40169380 0.0054568832082
186.290049920 0.0406384092110
42.2511763460 0.1802559388800
11.6765579320 0.4631512175500
3.59305064820 0.4408717331400
0 0 1 2 1
0.46866670000 1.0000000000000
0 0 1 0 1
0.23547146000 1.0000000000000
0 2 3 2 1
9.46809706210 0.0383878717280
2.01035451420 0.2111702511200
0.54771004707 0.5132817211400
0 2 1 0 1
0.21046214000 1.0000000000000
0 3 1 0 1
0.52857792000 1.0000000000000
8 5
0 0 6 2.0 1.0
0.5484671660D+04 0.1831074430D-02
0.8252349460D+03 0.1395017220D-01
0.1880469580D+03 0.6844507810D-01
0.5296450000D+02 0.2327143360D+00
0.1689757040D+02 0.4701928980D+00
0.5799635340D+01 0.3585208530D+00
0 1 3 6.0 1.0
0.1553961625D+02 -0.1107775495D+00 0.7087426823D-01
0.3599933586D+01 -0.1480262627D+00 0.3397528391D+00
0.1013761750D+01 0.1130767015D+01 0.7271585773D+00
0 1 1 0.0 1.0
0.2700058226D+00 0.1000000000D+01 0.1000000000D+01
0 1 1 0.0 1.0
0.8450000000D-01 0.1000000000D+01 0.1000000000D+01
0 3 1 0.0 1.0
0.8000000000D+00 0.1000000000D+01
1 4
0 0 3 1.0 1.0
0.1873113696D+02 0.3349460434D-01
0.2825394365D+01 0.2347269535D+00
0.6401216923D+00 0.8137573261D+00
0 0 1 0.0 1.0
0.1612777588D+00 0.1000000000D+01
0 2 1 0.0 1.0
0.2200000000D+01 0.1000000000D+01
0 2 1 0.0 1.0
0.5500000000D+00 0.1000000000D+01
99 0
END
DFT
EXCHANGE
PBE
CORRELAT
PBE
CHUNKS
200
XXLGRID
END
TOLDEE
9
TOLINTEG
11 11 11 11 20
SHRINK
20 20
END
With the optimized cell, perform the vibrational analysis
graphene
SLAB
1
2.47398205 2.47398694 120.000076
2
6 3.551099899926E-07 1.528176553945E-07 1.149651240873E-06
6 3.333329778900E-01 -3.333334868177E-01 -1.149651240873E-06
FREQCALC
NUMDERIV
2
STEPSIZE
0.003
END
END
6 6
0 0 5 2 1
1238.40169380 0.0054568832082
186.290049920 0.0406384092110
42.2511763460 0.1802559388800
11.6765579320 0.4631512175500
3.59305064820 0.4408717331400
0 0 1 2 1
0.46866670000 1.0000000000000
0 0 1 0 1
0.23547146000 1.0000000000000
0 2 3 2 1
9.46809706210 0.0383878717280
2.01035451420 0.2111702511200
0.54771004707 0.5132817211400
0 2 1 0 1
0.21046214000 1.0000000000000
0 3 1 0 1
0.52857792000 1.0000000000000
8 5
0 0 6 2.0 1.0
0.5484671660D+04 0.1831074430D-02
0.8252349460D+03 0.1395017220D-01
0.1880469580D+03 0.6844507810D-01
0.5296450000D+02 0.2327143360D+00
0.1689757040D+02 0.4701928980D+00
0.5799635340D+01 0.3585208530D+00
0 1 3 6.0 1.0
0.1553961625D+02 -0.1107775495D+00 0.7087426823D-01
0.3599933586D+01 -0.1480262627D+00 0.3397528391D+00
0.1013761750D+01 0.1130767015D+01 0.7271585773D+00
0 1 1 0.0 1.0
0.2700058226D+00 0.1000000000D+01 0.1000000000D+01
0 1 1 0.0 1.0
0.8450000000D-01 0.1000000000D+01 0.1000000000D+01
0 3 1 0.0 1.0
0.8000000000D+00 0.1000000000D+01
1 4
0 0 3 1.0 1.0
0.1873113696D+02 0.3349460434D-01
0.2825394365D+01 0.2347269535D+00
0.6401216923D+00 0.8137573261D+00
0 0 1 0.0 1.0
0.1612777588D+00 0.1000000000D+01
0 2 1 0.0 1.0
0.2200000000D+01 0.1000000000D+01
0 2 1 0.0 1.0
0.5500000000D+00 0.1000000000D+01
99 0
END
DFT
EXCHANGE
PBE
CORRELAT
PBE
CHUNKS
200
XXLGRID
END
TOLDEE
9
TOLINTEG
11 11 11 11 20
SHRINK
20 20
END
NUMDERIV set to 2 means double-sided finite difference while calculating the Hessian. STEPSIZE sets the step size for finite displacements in Å.
The calculation produces the main ouput file xxxx.out and the Hessian file HESSFREQ.DAT.
Then follow the instructions in quickstart example to do local mode analysis with LModeA-nano.