Relaxation of Structure

This tutorial will guide the user through the general workflow to perform a relaxation of a water dimer by using python based packages such as eDFTpy and the Atomic Simulation Environment (ASE).

Note

To avoid mixing some files, we strongly recommend setup each calculation (e.g. Optimization of density, relaxation of a structure of an MD simulation) in separate folders.

Input files

To carry out the relaxation of a structure, you have a set of input files divided in four different types:

• Global system coordinates

• eDFTpy input file

• Aditional subsystem configuration

• ASE input file

Additional to these input files and the same that for the previous SCF tutorials, you also need the pseudopotential files for each of the atomic species described in your global system coordinates input file.

Global System Coordinates

In this tutorial, we are going to use the same global system structure provided in the second optimization of densitytutorial composed of a binary salt (NaCl) molecule surrounded by a water dimer.

salt_h2o_2.xyz

8
Lattice="12.0 0.0 0.0 0.0 12.0 0.0 0.0 0.0 12.0" Properties=species:S:1:pos:R:3 pbc="T T T"
Na       2.00000000       2.50000000       4.2000000
Cl       6.00000000       7.50000000       9.20
O        2.00000000       2.00000000       2.11926200
H        2.00000000       2.76323900       1.52295300
H        2.00000000       1.23676100       1.52295300
O        6.00000000       6.00000000       6.11926200
H        6.00000000       6.76323900       5.52295300
H        6.00000000       5.23676100       5.52295300

eDFTpy input file

This input file follows the same structure that in the SCF Calculation with heterogeneous subsystems tutorial with some additional keywords in the Chloride subsystem section with the goal of setup the use of gaussian smearing to describe the occupations explained below:

input.ini

[JOB]
task               = Optdensity

[PATH]
pp                 = ./
cell               = ./

[PP]
O                  = o_pbe_v1.2.uspp.F.UPF
H                  = h_pbe_v1.4.uspp.F.UPF
Na                 = na_pbe_v1.5.uspp.F.UPF
Cl                 = cl_pbe_v1.4.uspp.F.UPF 

[OPT]
maxiter            = 200
econv              = 1e-6

[MOL]
charge-Na  = 1
charge-H2O = 0
charge-Cl = -1

[GSYSTEM]
cell-file          = salt_h2o_2.xyz
grid-ecut          = 1200
exc-x_str          = gga_x_pbe
exc-c_str          = gga_c_pbe
kedf-kedf          = GGA
kedf-k_str         = revAPBEK

[SUB_NA]
calculator         = qe
embed              = KE XC
cell-split         = 0.5 0.5 0.5
cell-index         = 0:1
grid-ecut          = 2400
mix-coef            = 0.7

[SUB_CL]
calculator         = qe
embed              = KE XC
cell-split         = 0.5 0.5 0.5
cell-index         = 1:2
basefile           = qe_in_cl.in
density-use_gaussians = True
mix-coef           = 0.7

[SUB_H2O]
calculator         = qe
embed              = KE XC
cell-split         = 0.5 0.5 0.5
cell-index         = 2:
decompose-method   = distance
decompose-radius-O = 0.90
decompose-radius-H = 0.60
grid-ecut          = 2400
mix-coef            = 0.7

• [SUB_CL]

basefile           = qe_in_cl.in           # File that contains additional QE keywords for Chloride subsystem
density-use_gaussians = True               # Using gaussian smearing

Aditional subsystem configuration

We defined a base file for the Chloride subsystem in the eDFTpy input file, we should also provide this input file with a Quantum Espresso input format. The goal of this input file is to provide an additional calculation setup for the determined subsystem such as LDA+U (see the additional QEpy configuration to enable this type of calculation), which means you can use one base file for each subsystem if you require and specify it in your eDFTpy input file.

qe_in_cl.in

&SYSTEM
   ecutwfc          = 40
   ecutrho          = 400
   nosym            = .true.
   occupations      = 'smearing'
   degauss          = 0.001
   smearing         = 'gaussian'
/
&ELECTRONS
/

ASE input file

The goal of this input file is to configure the relaxation of a strcuture driven by ASE software, merging the subsystem configuration provide by the eDFTpy driver calculation. This input file follows the ASE input format for a relaxation calculation.

ase_relax.py

#!/usr/bin/env python3

import os
import numpy as np
from edftpy.api.api4ase import eDFTpyCalculator
from edftpy.config import read_conf
from edftpy.interface import conf2init
from edftpy.mpi import graphtopo, sprint, pmi
np.random.seed(8888)
############################## initial ##############################
conf = read_conf('./input.ini')
conf2init(conf, pmi.size > 0)
cell_file = conf["PATH"]["cell"] +os.sep+ conf['GSYSTEM']["cell"]["file"]
#-----------------------------------------------------------------------

import ase.io
from ase.io.trajectory import Trajectory
from ase import units
from ase.md.verlet import VelocityVerlet
from ase.constraints import FixAtoms
from ase.md.andersen import Andersen
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.constraints import FixBondLength
from ase.optimize import BFGS, LBFGS, FIRE
from ase.optimize.sciopt import SciPyFminBFGS, SciPyFminCG

atoms = ase.io.read(cell_file)

calc = eDFTpyCalculator(config = conf, graphtopo = graphtopo)
atoms.set_calculator(calc)

trajfile = 'opt.traj'
af = atoms
opt = LBFGS(af, trajectory = trajfile, memory = 10, use_line_search = False)
# opt = SciPyFminCG(af, trajectory = trajfile)

opt.run(fmax = 0.3)

traj = Trajectory(trajfile)
ase.io.write('opt.xyz', traj[-1])

To setup your own calculation you should take into account and reconfigure the next keywords sections:

conf = read_conf('./input.ini')                                                 # Change if your eDFTpy input file has a different name

opt = LBFGS(af, trajectory = trajfile, memory = 10, use_line_search = False)    # In this tutorial we use the limited memory version of the BFGS algorithm, check ASE wwebsite to use a different algorithm

opt.run(fmax = 0.3)

Note

In this tutorial we use the limited memory version of the BFGS algorithm, check ASE website to use a different algorithm. On the other hand. we use a big value for fmax, you should change for a real calculation setup.

Running eDFTpy

Once you finish setting out the calculation, you can run the calculation using:

$ mpirun -n 4 python -m ase_relax.py > log

Note

Here the log file storages the summary of the calculation.

Warning

If you are using python3 version you should change all the python command lines wrote here from python to python3

Output files

Once your calculation finishes you will have the same set of output files from eDFTpy calculations with additional files as a product of the relaxation of the global system.

Additional output Files

• edftpy_gsystem.xyz Global system coordinates for the last step of the simulation

• opt.xyz Coordinates for the global system optimized structure

• opt.traj Trajectory file for the relaxation proccess

Preliminary Analysis

Once you finished your simulation you would like to analyze your results, most of the time your structural analysis will require your trajectory file (opt.traj). Although this trajectory formatis only useful if you are using the ASE package to perform visual analysis.

Open the trajectory file with the ASE visualization tool by typing:

$ ase gui opt.traj

Once you open the relaxation trajectory file using the visualization tool included in ASE package you will observe the step by step structure change and a graphic generator which includes a set of graphics check ASE’s GUI. Here we show the variation of the total energy (eV) in function of the number steps for two different relaxations of the binary salt (NaCl) molecule surrounded by a water dimer, using the fmax = 0.3 value provided in this tutorial (left) and a fmax = 0.05 (Right).

Total energy Variation

Convert the format of a trajectory file

If you would like to use any other package to carry out structural analysis will be useful for you to covert this trajectory file to a most common type as xyz. We developed a python script that you can use to convert your trajectory files by typing:

$ python -m edftpy –-convert –-traj opt.traj -o opt-traj.xyz

You can check the command line miscellaneous for eDFTpy in the command line section of the sofware website.

Restarting a Relaxation of a structure

Restarting a calculation could be a powerful option taking into account the time-consuming disadvantage of an ab-initio molecular dynamics simulation. Here we give the user a guide on how can restart a relaxation that for a pecific reason did not finish or should be elongated.

To restart a relaxation calculation the user should modify the eDFTpy and ASE input files either Structure optimizartion or MD production, we suggest the next six-step workflow to have succes restarting a calculation.

1. Copy your trajectory file (opt.traj) to a new file:

   $ cp opt.traj  opt-1.traj
   

2. Create a new global system coordinate file using the last geometry coordinates:

   $ cp edftpy_gsystem.xyz rst-1.xyz
   

3. Change the file name of the cell-file in the eDFTpy input file

   • from

   cell-file          = h2o_2.xyz
   

   • to

   cell-file          = rst-1.xyz
   

4. Change the cell file in the ASE input file to your old trajectory file

   • from

   atoms = ase.io.read(cell_file)
   

   • to

   atoms = ase.io.read("opt-1.traj")
   

5. Restart your simulation using:

   $ mpirun -n 4 python -m ase_relax.py > log
   

6. Once you fininsh your structure optimization after one of many restartings you could concatenate the trajectory files using:

   $ python -m edftpy --convert --traj opt-1.traj opt-2.traj -o opt-tot.traj