Structure of the water molecule¶
This Notebook is oriented on Chapter 2 of the CP-PAW Hands-On Course on First-Principles Calculations. It covers the same content while showing different ways how to use the cppawAse interface.
cppawAse preparation¶
First the path to the cppawAse Files has to be set and the main module loaded. Within the jupyter notebook a path to the directory for the calculation has to be speciefied and added to the calculation label.
In [1]:
import sys
acPath = '/home/pmk/paw/hiwi/AseCppaw'
sys.path.append(acPath)
import AseCppaw
import AseCppawInputFiles
path0 = acPath + '/tests/h2o/'
Input¶
Within the ASE simulations are run by a calculator that is set to an atoms object containing atom coordinates, cell and border conditions. In cppawAse the calculator is the cppaw class it contains the inputfiles (at least one structure and one control file) and parameters how to run the calculation (single-core, parallel on a cluster with queueing system), backups, restart and what to log. It also contains the results of the simulation. Like total energy, band gap, fermi level etc. The instructions list can minimally contain one entry consisting of the keys 'strc' und 'cntl' or many entries that are simulated in order of the list.
The CP-PAW way¶
In [2]:
# loading the input files
cntl = AseCppawInputFiles.loadInputFile(path0 + 'h2o.cntl1')
strc = AseCppawInputFiles.loadInputFile(path0 + 'h2o.strc')
# setting them as cppaw input, label is the subfolder were the simulation is located
calc = AseCppaw.cppaw(instructions=[{'cntl': cntl, 'strc': strc}],
label=path0 + 'h2o_220530',
report='l',
logCalc='l')
In [3]:
# the ase atoms object is automatically created from the structure file and can be viewed
from ase.visualize import view, ngl
v = view(strc.atoms, viewer='x3d')
v
Out[3]:
The ASE way¶
In [4]:
# the ASE provides some common molecules as atoms object within the molecule module
from ase.build import molecule
h2o = molecule('H2O')
v = view(h2o, viewer='x3d')
v
Out[4]:
In [5]:
# contrary to the CP-PAW structure the h2o molecule doesn't have a unit cell
# the recomended 6.0 Angstrom distance between periodic images can be archived with ASE atom functions
h2o.center(vacuum=6.0)
# setting the periodic boundary condition pbc to False while trigger the CP-PAW !ISOLATE branch
h2o.set_pbc(False)
In [6]:
# ASE automatically chooses an bcc lattice that is suboptimal for CP-PAW calculations
# the cell can be set manually as an list of lists or numpy array
h2o = molecule('H2O')
h2o.set_cell([[0.0, 6.0, 6.0], [6.0, 0.0, 6.0], [6.0, 6.0, 0.0]])
# without viewer='x3d' defined the structure including the unit cell is shown in an interactive external window
view(h2o)
Out[6]:
In [8]:
# for more complicated cases the lattice parameters can be determined automatically for a given celltype and minimal
# distance between periodic images. A structure file has to be already present.
import AseCppawTools
h2o = molecule('H2O')
cellparameter = AseCppawTools.makeLattice(celltype='FACE CENTERED', dist=7.0, filename='h2o', directory=path0,
flags='-c')
h2o.set_cell(cellparameter['cell'])
view(h2o)
Out[8]:
In [9]:
# the calculater can now be initialized with just the atoms object as input using default control files and setups
calcAse = AseCppaw.cppaw(atoms=h2o,
label=path0 + 'h2o_220915',
report='l',
logCalc='l')
interacting with input objects¶
In [10]:
# parameters can be accessed by entering the branch as list ...
T = strc.getValue(['!LATTICE', 'T'])
print(T)
In [11]:
# The automatically generated input files can be accesed in the instructions in the parameters
# of the ase calculator object
aseStrc = calcAse.parameters['instructions'][0]['strc']
In [12]:
# ... or as string seperated by , or ; or blank space list
T = aseStrc.getValue('!LATTICE;T')
print(T)
In [13]:
# parameters can be changed in a simular way
print(cntl.getValue('!GENERIC;DT'))
cntl.setValue('!GENERIC;DT', 5.)
print(cntl.getValue('!GENERIC;DT'))
In [14]:
# if multiple branches exist a list of lists is used where the first list is the parent branch
# the second the leaf value pair to identify the branch and the last list the branches from the parent branch
# the identifier must be a list of the unique leaf here 'NAME' and it's value
Oid = strc.getValue([['!SPECIES'], ['NAME','O_'], ['ID']])
Oid
Out[14]:
In [15]:
# alternatively the lists can be indicated in a string with pipe | characters
Hid = strc.getValue('!SPECIES|NAME;H_|ID')
Hid
Out[15]:
In [16]:
# if usage or meaning of a parameter is unclear get_help returns the description in the CP-PAW manual
from IPython.display import display, Latex
display(Latex(strc.getHelp('!LATTICE;T')))
Running a simulation¶
In [17]:
# the calculation can be started by calling any function that returns a property that requires a calculation
Eurlx = calcAse.get_potential_energy()
# The energy is given in electron volts, the standard unit of the ASE
print(Eurlx)
In [18]:
# the most recent results are stored in the results dict
results = calcAse.results
for key in results.keys():
print(key, results[key])
In [19]:
# the progress can be plotted while using for opt the same letters as in the paw_show tool
calcAse.protocol.showGraph(opt='ec')
In [20]:
# it can also display in kj/mol or ev and only track the last n steps
calcAse.protocol.showGraph(opt='ec', unit='kj/mol', nSteps=50)
In [21]:
# to run the structure relaxation we change the parameters of the default control files
calcAse.parameters
Out[21]:
In [24]:
# we remove the default start control file
cntlDict = calcAse.parameters['instructions'].pop(1)
cntlDict = calcAse.parameters['instructions'].pop(1)
cntlDict
In [28]:
# and change the electron relaxation control file to include rdyn
cntlStrg = cntlDict['cntl']
rlxCntl = AseCppawInputFiles.inputtoobject(cntlStrg)
rlxCntl.setValue(['!RDYN'], True)
rlxCntl.setValue(['!PSIDYN', '!AUTO', 'FRIC(-)'], 0.0)
calcAse.parameters['instructions'].append({'cntl': rlxCntl.getInputfileStr()})
In [29]:
# and start the calculation again with the changed parameters
calcAseRlxr = AseCppaw.cppaw(label=path0 + 'h2o_220915', **calcAse.parameters)
calcAseRlxr.calculate()
In [30]:
# all results of the calculation can be directly accesed witin the result attribute
print(calcAseRlxr.results)
In [ ]: