[ase-users] Thermochemistry with cp2k and ASE
Hasan Al-Mahayni
hasanalmahayni at gmail.com
Sun Jun 14 18:50:45 CEST 2020
Hi all,
I am trying to run a gibbs free energy calculation (vibrational analysis)
on an optimized copper slab with CO adsorbed on it. I get the
following error:
ValueError: Imaginary vibrational energies are present.
I attached the python code(harmonic) and the xyz slab (optslab).
Thanks,
Hasan.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20200614/b1afcf4d/attachment.html>
-------------- next part --------------
50
C 2.5524988134 4.4240144500 6.1878260399
O 2.5596000724 4.4387853693 7.3447274640
Cu -0.0000510000 1.4757220000 0.0000000000
Cu -1.2780390000 3.6892620000 0.0000000000
Cu -2.5560260000 5.9028020000 0.0000000000
Cu -3.8340140000 8.1163410000 0.0000000000
Cu 2.5559240000 1.4757220000 0.0000000000
Cu 1.2779360000 3.6892620000 0.0000000000
Cu -0.0000510000 5.9028020000 0.0000000000
Cu -1.2780390000 8.1163410000 0.0000000000
Cu 5.1118980000 1.4757220000 0.0000000000
Cu 3.8339110000 3.6892620000 0.0000000000
Cu 2.5559230000 5.9028020000 0.0000000000
Cu 1.2779360000 8.1163410000 0.0000000000
Cu 7.6678730000 1.4757220000 0.0000000000
Cu 6.3898860000 3.6892620000 0.0000000000
Cu 5.1118980000 5.9028020000 0.0000000000
Cu 3.8339110000 8.1163410000 0.0000000000
Cu 1.2701433807 0.7197902303 2.1103978801
Cu -0.0104990319 2.9331876305 2.1144558486
Cu -1.2869625107 5.1465145876 2.1099729454
Cu -2.5648702151 7.3594544757 2.1112349011
Cu 3.8224161241 0.7197508490 2.1096256991
Cu 2.5479069140 2.9490270924 2.1320858162
Cu 1.2821226513 5.1400908674 2.1325640937
Cu -0.0111458421 7.3576780300 2.1101718226
Cu 6.3796861662 0.7194193582 2.1122926294
Cu 5.1029660263 2.9337761965 2.1138200524
Cu 3.8146274171 5.1402937707 2.1323172520
Cu 2.5471354576 7.3584692997 2.1140797306
Cu 8.9378104585 0.7190902209 2.1133669911
Cu 7.6581722369 2.9342219530 2.1126240941
Cu 6.3803804998 5.1471625374 2.1100607339
Cu 5.1048430765 7.3581023926 2.1110490069
Cu -0.0142149943 -0.0282370127 4.2238112902
Cu -1.2944705363 2.1842796929 4.2224015505
Cu -2.5698870924 4.3981381667 4.2216335926
Cu -3.8467390127 6.6119196259 4.2219324670
Cu 2.5427846843 -0.0302658152 4.2232017343
Cu 1.2497913885 2.1661295371 4.2148110446
Cu -0.0377630003 4.3965961400 4.2150970929
Cu -1.2941376227 6.6131419190 4.2230901261
Cu 5.0979895051 -0.0283877994 4.2234595804
Cu 3.8350698300 2.1646613493 4.2162219857
Cu 2.5405208126 4.3964551900 4.3232403598
Cu 1.2521343604 6.6348495231 4.2162849702
Cu 7.6528165086 -0.0250634030 4.2241134653
Cu 6.3795201840 2.1826969567 4.2236706783
Cu 5.1221917665 4.3958834536 4.2152577526
Cu 3.8307408306 6.6338514313 4.2149119860
-------------- next part --------------
#imports
import numpy as np
from ase.build import molecule
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo
from ase.calculators.cp2k import CP2K
from ase.io import read
from ase import io
inp= """
! &EXT_RESTART
! RESTART_FILE_NAME cp2k-1.restart
! &END
&MOTION
&PRINT
&RESTART
&EACH
&END EACH
&END RESTART
&TRAJECTORY
FORMAT XYZ
&END TRAJECTORY
&CELL
&END CELL
&END PRINT
&END MOTION
&FORCE_EVAL
&DFT
UKS .TRUE.
CHARGE 0
MULTIPLICITY 1
WFN_RESTART_FILE_NAME cp2k-RESTART.kp
&MGRID
NGRIDS 5
RELATIVE_CUTOFF 50
&END MGRID
&QS
EPS_DEFAULT 1.0E-12
METHOD GPW
EXTRAPOLATION USE_GUESS
&END QS
&SCF
SCF_GUESS RESTART
ADDED_MOS 200
CHOLESKY OFF
&SMEAR T
METHOD FERMI_DIRAC
ELECTRONIC_TEMPERATURE 3.0000E+2
&END SMEAR
&MIXING T
METHOD BROYDEN_MIXING
ALPHA 4.0000000000000002E-01
BETA 1.5000000000000000E+00
NMIXING 5
NBUFFER 8
&END MIXING
&END SCF
&POISSON
POISSON_SOLVER PERIODIC
PERIODIC XYZ
&END POISSON
&KPOINTS
SCHEME MONKHORST-PACK 8 8 8
FULL_GRID .TRUE.
&END KPOINTS
&END DFT
&END FORCE_EVAL
"""
slab=io.read('slab.xyz')
slab.center(vacuum=15.0)
slab.set_calculator(CP2K(xc='PBE',cutoff=500 ,max_scf=200,poisson_solver='none' ,inp=inp))
#potentialenergy = slab.get_potential_energy() #energy optimized
potentialenergy=-64440.33054394734
indices = [] #specify which are adsorbate atoms by index number
for i in range(0,2):
indices.append(i)
#this section calculates the vibrational energy
vib = Vibrations(slab,indices=indices,name='vib')
vib.run()
vib.summary(log='vib.txt')
for mode in range(len(indices)*3):
vib.write_mode(mode)
vib_energies = vib.get_energies()
#this section inputs the necessary input to calculate thermodynamic properties
thermo = HarmonicThermo(vib_energies=vib_energies, potentialenergy=potentialenergy )
G = thermo.get_gibbs_energy(temperature=298.15, pressure=101325.)
More information about the ase-users
mailing list