# ATST-Tools
Advanced ASE Transition State Tools for ABACUS, including:
- NEB, including CI-NEB, IT-NEB and others
- AutoNEB: an automatic NEB workflow
- Dimer
- Vibration analysis

Version v1.3.1

Copyright @ QuantumMisaka from PKU & AISI

## Dependencies:
- [ASE](https://wiki.fysik.dtu.dk/ase/about.html)
- [ABACUS](https://abacus.deepmodeling.com/en/latest/)
- [ASE-ABACUS interface](https://gitlab.com/1041176461/ase-abacus)
- [GPAW](https://wiki.fysik.dtu.dk/gpaw/install.html) if one wants to run NEB images relaxation in parallel

Notice: GPAW and ABACUS should be dependent on same MPI and libraries environments. 
For instance, if your ABACUS is installed by Intel-OneAPI toolchain, your GPAW should NOT be dependent on gcc-toolchain like OpenMPI and OpenBLAS.

ATST-Tools is Under actively development, please let me know if any problem occurs.

## Workflow

![ATST-workflow](./img/ATST-workflow.png)

## Workflow libraries and setting
All workflow library files and re-constructed ASE libraries will be put in `./source` directory. including:
```bash
source
├── abacus_autoneb.py
├── abacus_dimer.py
├── abacus_neb.py
└── my_autoneb.py
```
Before use running scripts, you should add these libraries into your PYTHONPATH:
```bash
export PYTHONPATH=/path/to/source:$PYTHONPATH
```

## Developing
- [x] Use interface to read ABACUS STRU file and ABACUS output
- [x] Flexible input for different NEB method in ASE
- [x] `DyNEB` implementation and test
- [x] Now used optimum option: idpp-guess + IT-NEB + CI-NEB parallel method
- [x] Give an initial guess print-out of NEB images
- [x] Decoupling to init-guess -> NEB calculation -> result post-process
- [x] Parallel computing during images relaxation by `gpaw python`
- [x] `AutoNEB` implementation
- [x] Give user a way to specify constraint information (along one direction) and initial magnetic moments (for each atom of one element) for NEB initlal guess chain and two end points (while these information cannot be read from calculation result)
- [x] More test in surface reaction system
- [x] Connected to Dimer method
- [x] More test in magnetic surface reaction system
- [x] Put calculation setting in an independent file (decoupling *_run.py)
- [x] Vibration analysis scripts 
- [x] More easily and proper way to do continuation calculation for NEB and AutoNEB. (next target)
- [ ] More fiexible options for NEB, Dimer and AutoNEB, like full properties in trajectory file, and fiexibly utilize SCF wavefunction/charge output files from previous calculation.


## NEB workflow

### Method
- For serial NEB calculation, DyNEB, namely dynamic NEB method `ase.mep.neb.DyNEB` is for default used.
- For parallel NEB calculation, `ase.mep.neb.NEB` traditional method is for default used.
- The Improved Tangent NEB method `IT-NEB` and Climbing Image NEB method `CI-NEB` in ASE are also default used in this workflow, which is highly recommended by Sobervea. In `AutoNEB`, `eb` method is used for default, but Improved Tangent method is also recommended.
- Users can change lots of parameter for different NEB setting. one can refer to [ASE NEB calculator](https://wiki.fysik.dtu.dk/ase/ase/neb.html#module-ase.neb) for more details: 
- The workflow also support use `AutoNEB` method in ASE。You can view AutoNEB method in paper below. Also. One can refer to [AutoNEB](https://wiki.fysik.dtu.dk/ase/ase/neb.html#autoneb) to view it. 
> E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys, 145, 094107, 2016. (doi: 10.1063/1.4961868)

The AutoNEB method in ASE lies in `ase.mep.autoneb.AutoNEB` object, which will do NEB calculation in following steps:
1. Define a set of images and name them sequentially. Must at least have a relaxed starting and ending image. User can supply intermediate guesses which do not need to have previously determined energies (probably from another
NEB calculation with a lower level of theory)
1. AutoNEB will first evaluate the user provided intermediate images
2. AutoNEB will then add additional images dynamically until n_max is reached
3. A climbing image will attempt to locate the saddle point
4. All the images between the highest point and the starting point are further relaxed to smooth the path
5. All the images between the highest point and the ending point are further relaxed to smooth the path

Step 4 and 5-6 are optional steps. Note that one can specify different `fmax` for CI-NEB in step 4 compared with other NEB calculation, which can be set as `fmax=[fmax1, fmax2]` in `AutoNEB` object.

> Notice: in surface calculation and hexangonal system, the vaccum and c-axis should be set along y-direction but not z-direction, which is much more efficient for ABACUS calculation.


### Usage
#### Basic NEB
The NEB workflow is based on 3 main python scripts and 1 workflow submit script. Namely:

- `neb_make.py` will make initial guess for NEB calculation, which is based on ABACUS (and other calculator) output files of initial and final state. This script will generate `init_neb_chain.traj` for neb calculation. Also, You can do continuation calculation by using this script. You can get more usage by `python neb_make.py`. 
- `neb_run.py` is the key running script of NEB, which will run NEB calculation based on `init_neb_chain.traj` generated by `neb_make.py`. This script will generate `neb.traj` for neb calculation. Users should edit this file to set parameters for NEB calculation. sereal running can be performed by `python neb_run.py`, while parallel running can be performed by `mpirun gpaw python neb_run.py`.
When running, the NEB trajectory will be output to `neb.traj`, and NEB images calculation will be doing in `NEB-rank{i}` directory for each rank which do calculation of each image. 
- `neb_post.py` will post-process the NEB calculation result, which will based on `neb.traj` from neb calculation. This script will generate nebplots.pdf to view neb calculation result, and also print out the energy barrier and reaction energy. You can get more usage by `python neb_post.py`. Meanwhile, users can also view result by `ase -T gui neb.traj` or `ase -T gui neb.traj@-{n_images}:` by using ASE-GUI
- `neb_submit.sh` will do all NEB process in one workflow scripts and running NEB calculation in parallel. Users should edit this file to set parameters for NEB calculation. Also this submit script can be used as a template for job submission in HPC. the Default setting is for `slurm` job submission system.

#### AutoNEB method 
In ATST-Tools, the AutoNEB method can be easily used by the following scripts
- `autoneb_run.py` is the key running script for `AutoNEB` method, which is like `neb_run.py` but the NEB workflow in `AutoNEB` is enhanced and the I/O logic have some difference. Users can use it with `mpirun gpaw python autoneb_run.py` by existing `init_neb_chain.traj` which can only contain initial and final state or contain some initial-guess.
- `autoneb_submit.sh` will do all AutoNEB process in one workflow and running AutoNEB calculation in parallel. Users should edit this file to set parameters for AutoNEB calculation. Also this submit script can be used as a template for job submission in HPC. the Default setting is for `slurm` job submission system.
- `neb_make.py` and `neb_post.py` can be used for `AutoNEB` method, but the workflow have slight difference. 

#### Running
Users can run NEB each step respectively: 
1. `python neb_make.py [INIT/result] [FINAL/result] [n_max]` to create initial guess of neb chain
   1. Also You can use `python neb_make.py -i [input_traj_file] [n_max]` to create initial guess from existing traj file, which can be used for continuation calculation.
2. `python neb_run.py` or `mpirun -np [nprocs] gpaw python neb_run.py` to run NEB calculation
3. `python neb_post.py neb.traj [n_max]` to post-process NEB calculation result

Users can run AutoNEB each step respectively:
1. `python neb_make.py [INIT/result] [FINAL/result] [nprocs]` to create initial guess of neb chain
2. `mpirun -np [nprocs] gpaw python autoneb_run.py` to run AutoNEB calculation
3. `python neb_post.py --autoneb run_autoneb???.traj` to post-process NEB calculation result

Also, user can run each step in one script `neb_submit.sh` by `bash neb_submit.sh` or `sbatch neb_submit.sh`. AutoNEB scripts usage is like that. 

> Notice: Before you start neb calculation process, make sure that you have check the nodes and cpus setting and other setting like n_max, constraints and initial magnetic moments in `*neb_submit.sh` and `*neb_run.py` to make sure that you can reach the highest performance and reach the simulation result you want !!!   

#### Visualize
You can simply use ASE-GUI to have a view of NEB or AutoNEB trajectory.

For NEB, you can view all running trajectory by 
```bash
ase -T gui neb.traj 
```
You can also view the last 10 images by
```bash
ase -T gui neb.traj@-10:
```
For AutoNEB, the most recent NEB path can always be monitored by:
```bash
ase -T gui -n -1 run_autoneb???.traj
```

#### Continuation calculation for NEB
If NEB or AutoNEB is break down somehow, you can do continuation calculation based on saved trajectory files and ATST-Tools scripts.

For NEB, you can simply:
```bash
python neb_make.py -i neb.traj [n_max] [fix and mag information]
```
to generate `init_neb_chain.traj` for continuation calculation. You can also `python neb_post.py neb.traj` to generate the latest neb band `neb_latest.traj` and do continuation calculation by `python neb_make.py -i neb_latest.traj [n_max]`. note that `n_max = n_image - 2`

For AutoNEB, you need to get `neb_latest.traj` in a more compicated way, especially for the first NEB step in AutoNEB running:
```bash
python traj_collect.py ./AutoNEB_iter/run_autoneb???iter[index].traj
```
to generate `collection.traj` from certain index (like 006) stage of AutoNEB calculation. Another way to do the same thing is:
```bash
python traj_collect.py ./run_autoneb???.traj
```
or
```bash
python traj_collect.py ./AutoNEB_iter/run_autoneb???iter00[i].traj
```
or
```bash
python traj_collect.py ./AutoNEBrun_rank?/STRU
```
When `collection.traj` is gotten, one can do
```bash
python neb_make.py -i collection.traj [n_max] [fix and mag information]
```
to generate `init_neb_chain.traj` for continuation calculation.

If one want to continue calculation with the interrupted autoneb, one can:
```bash
python traj_collect.py --no-calc ./run_autoneb???.traj 
```
and `--no-calc` noted can be anywhere.

> Note: Linux shell will automatically detect and sort number of index, so you will not be worried about using format like `run_autoneb???iter005.traj`, the consequence will be right, for example:
```bash
 test> ll run_autone*.traj
-rw-r--r-- 1 james james 6.0K Nov 24 20:35 run_autoneb000.traj
-rw-r--r-- 1 james james 6.4K Nov 24 20:35 run_autoneb001.traj
-rw-r--r-- 1 james james 6.4K Nov 24 20:35 run_autoneb002.traj
-rw-r--r-- 1 james james 531K Nov 24 20:35 run_autoneb003.traj
-rw-r--r-- 1 james james 531K Nov 24 20:35 run_autoneb004.traj
-rw-r--r-- 1 james james 531K Nov 24 20:35 run_autoneb005.traj
-rw-r--r-- 1 james james 531K Nov 24 20:35 run_autoneb006.traj
-rw-r--r-- 1 james james 6.4K Nov 24 20:35 run_autoneb007.traj
-rw-r--r-- 1 james james 6.4K Nov 24 20:35 run_autoneb008.traj
-rw-r--r-- 1 james james 6.5K Nov 24 20:35 run_autoneb009.traj
-rw-r--r-- 1 james james 6.5K Nov 24 20:35 run_autoneb010.traj
-rw-r--r-- 1 james james 6.5K Nov 24 20:35 run_autoneb011.traj
-rw-r--r-- 1 james james 6.5K Nov 24 20:35 run_autoneb012.traj
-rw-r--r-- 1 james james 531K Nov 24 20:37 run_autoneb025.traj
```


#### Other scripts
Because ATST is originally based on ASE, the trajectory file can be directly read, view and analysis by `ase gui` and other ASE tools. Abide by `neb_make.py` and `neb_post.py`, We also offer some scripts to help you:
- `neb_dist.py`: This script will give distance between initial and final state, which is good for you to check whether the atoms in two image is correspondent, and is also a reference for setting number of n_max
- `traj_transform.py`: This script can transfer traj files into other format like `extxyz`, `abacus`(STRU), `cif` and so on (coming soon). Also if user specify `--neb` option, this script will automatically detect and cut the NEB trajectory when doing format transform. This script will be helpful for analysis and visualization of NEB trajectory.
- `traj_collect.py`: This script can collect structure files into a trajectory file, which is specifically used for NEB continuation calculation.


## Dimer workflow

### Method
(Waiting for update)


### Usage
The Dimer workflow is based on 2 main python scripts and 2 workflow submit script. Namely:
- `neb2dimer.py` can be used by `python neb2dimer [neb.traj] ([n_max])`, which will transform NEB trajetory `neb.traj` or NEB result trajectory `neb_result.traj` to Dimer input files,  including:
- - `dimer_init.traj` for initial state of Dimer calculation, which is the highest energy image, namely, TS state. 
- - `displacement_vector.npy` for displacement vector of Dimer calculation, which will be generated from position minus of the nearest image before and after TS point, and be normalized to 0.01 Angstrom. 
- `dimer_run.py` is the key running script of Dimer calculation, which will run Dimer calculation based on `dimer_init.traj` and `displacement_vector.npy` generated by `neb2dimer.py` or based on other setting. This script will generate `dimer.traj` for Dimer calculation trajectory. Users should edit this file to set parameters for Dimer calculation, and run Dimer calculation by `python dimer_run.py`. When running, any Dimer images calculation will be doing in `Dimer` directory.
- `dimer_submit.sh` will do Dimer workflow in one scripts. The Default setting is for `slurm` job submission system.

> Notice: this dimer method will have some problem in running process.


## Vibration Analysis 
The vibration analysis is based on `ase.vibrations.Vibrations` object, which can be used by `python vib_analysis.py` to do vibration analysis by finite displacement method for initial, final and transition state. The result will be printed out  and saved in `running_vib.out` file. All force matrix for displaced and normal mode will also be saved and printed.

Also, thermodynamic analysis will be performed based on `ase.thermochemistry.HarmonicThermo` object based on vibration analysis result and specified temperature.

## Relaxation
Relaxation method offer by ASE can be used by scripts from `relax` directory, which use ABACUS as SCF calculator. 


## Notices
### Property Loss in Trajectory 
Some property should be get via specific way from trajectory files, and some will be lost in trajetory files, 
- Stress property will not be stored in trajetory file
- In NEB calculation, the Force property for fixed atoms and Stress property will NOT be stored in trajectroy file, which is still in work
- in Dimer calculation, the Energy, Forces and Stress property will be stored in trajetory file after specified which properties need to be stored.
- in AutoNEB calculation, all property in processing trajectory will be stored in AutoNEB_iter directory, but in the result `run_autoneb???.traj`, the forces and stress information will be lost.


## Examples
- Li-diffu-Si: Li diffusion in Si, very easy example for serial and parallel NEB calculation
- H2-Au111: H2 dissociation on Au(111) surface. which will have NEB, AutoNEB and Dimer example. The barrier is around 1.1 eV consistent with existing paper and calculation result.
- N2-Cu111 : N2 dissociation on Cu(111) surface. which have high barrier as 4.0 eV but it's not hard to calculate
- CO-Pt111 : CO dissociation on Pt(111) surface. which have high barrier as 3.6 eV and including diffusion part. AutoNEB will always fail due to the low starting image. 
- Cy-Pt_graphene: Cyclohexane dehydrogenation on Pt-doped graphene surface. The barrier is around 1.3 eV. Noted that the `IT-NEB` result is wrong, but which is consistent to the result in VTST-Tools when using 4 image to do IT-NEB calculation.
- C-C coupling on Fe5C2(510) surface: C-C coupling on Fe5C2(510) surface with CO co-adsorption, which have 0.48 eV barrier in calculation result and existing paper. This is a typical example for reaction on magnetic surfaces.

More examples is welcomed from users. 

# JamesMisaka in 2023-11-27
# Run NEB calculation on NEB images by ASE-ABACUS
# part of ATST-Tools scripts

from ase.optimize import FIRE, BFGS
from ase.io import read, write
from abacus_neb import AbacusNEB

# setting
mpi = 8
omp = 4
fmax = 0.05  # eV / Ang
neb_optimizer = FIRE # suited for CI-NEB
neb_directory = "NEBrun"
algorism = "improvedtangent" # IT-NEB is recommended
climb = True
dyneb = False  
parallel = True
k = 0.10 # eV/Ang^2, spring constant
init_chain = "init_neb_chain.traj"
abacus = "abacus"
#lib_dir = "/lustre/home/2201110432/example/abacus"
lib_dir = ""
pseudo_dir = f"{lib_dir}"
basis_dir = f"{lib_dir}"
pp = {
      'C':'C_ONCV_PBE-1.0.upf',
      'H':'H_ONCV_PBE-1.0.upf',
      'Pt':'Pt_ONCV_PBE-1.0.upf',
      }
basis = {
         'C': 'C_gga_7au_100Ry_2s2p1d.orb',
         'H': 'H_gga_6au_100Ry_2s1p.orb',
         'Pt': 'Pt_gga_7au_100Ry_4s2p2d1f.orb'
         ,}
kpts = [2, 1, 2]
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'dft_functional': 'pbe',
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'd3_bj',
    'smearing_method': 'gaussian',
    'smearing_sigma': 0.001,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'scf_thr': 1e-6,
    'scf_nmax': 100,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'init_wfc': 'atomic',
    'init_chg': 'atomic',
    'cal_force': 1,
    'cal_stress': 1,
    'out_stru': 1,
    'out_chg': 0,
    'out_mul': 0,
    'out_wfc_lcao': 0,
    'out_bandgap': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
    'efield_pos_max': 0.0
}


if __name__ == "__main__":
# running process
# read initial guessed neb chain
    init_chain = read(init_chain, index=':')
    neb = AbacusNEB(init_chain, parameters=parameters, parallel=parallel,
                    directory=neb_directory, mpi=mpi, omp=omp, abacus=abacus, 
                    algorism=algorism, k=k, dyneb=dyneb)
    neb.run(optimizer=neb_optimizer, climb=climb, fmax=fmax)

# JamesMisaka in 2023-11-27
# Run AutoNEB calculation by ASE-ABACUS
# part of ATST-Tools scripts


from ase.optimize import FIRE, BFGS
from ase.io import read, write
from ase.parallel import world, parprint, paropen
#from pathlib import Path
from abacus_autoneb import AbacusAutoNEB

# setting
mpi = 16
omp = 4
neb_optimizer = FIRE # suited for CI-NEB
neb_directory = "AutoNEBrun"
# algorism = 'eb' # default for AutoNEB
algorism = "improvedtangent" # IT-NEB
init_chain = "init_neb_chain.traj"
climb = True
fmax = [0.20, 0.05]  # eV / Ang, 2 fmax for others and last CI-NEB in all-images
n_simul = world.size # only for autoneb, number of simultaneous calculation
n_images = 10 # only for autoneb, max number of all image, which should be reached
smooth_curve = False # True to do more neb step to smooth the curve.
k = 0.10 # eV/Ang^2, force constant of spring, 0.05 is from VTST-Tools
abacus = "abacus"
#lib_dir = "/lustre/home/2201110432/example/abacus"
lib_dir = ""
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
      'C':'C_ONCV_PBE-1.0.upf',
      'H':'H_ONCV_PBE-1.0.upf',
      'Pt':'Pt_ONCV_PBE-1.0.upf',
      }
basis = {
         'C': 'C_gga_7au_100Ry_2s2p1d.orb',
         'H': 'H_gga_6au_100Ry_2s1p.orb',
         'Pt': 'Pt_gga_7au_100Ry_4s2p2d1f.orb'
         ,}
kpts = [2, 1, 2]
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'dft_functional': 'pbe',
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'd3_bj',
    'smearing_method': 'gaussian',
    'smearing_sigma': 0.001,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'scf_thr': 1e-6,
    'scf_nmax': 100,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'cal_force': 1,
    'cal_stress': 1,
    'init_wfc': 'atomic',
    'init_chg': 'atomic',
    'out_stru': 1,
    'out_chg': 0,
    'out_mul': 0,
    'out_wfc_lcao': 0,
    'out_bandgap': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
    'efield_pos_max': 0.0
}


if __name__ == "__main__": 
# running process
# read initial guessed neb chain
    init_chain = read(init_chain, index=':')
    neb = AbacusAutoNEB(init_chain, parameters, algorism=algorism, 
                        directory=neb_directory, k=k,
                        n_simul=n_simul, n_max=n_images, 
                        abacus=abacus,  mpi=mpi, omp=omp, )
    neb.run(optimizer=neb_optimizer, climb=climb, 
                fmax=fmax, smooth_curve=smooth_curve)

Loading mkl version 2023.0.0
Loading tbb version 2021.8.0
Loading compiler-rt version 2023.0.0
Loading mpi version 2021.8.0
Loading compiler version 2023.0.0
Loading oclfpga version 2023.0.0
  Load "debugger" to debug DPC++ applications with the gdb-oneapi debugger.
  Load "dpl" for additional DPC++ APIs: https://github.com/oneapi-src/oneDPL

Loading abacus/3.4.2-icx
  Loading requirement: tbb/latest compiler-rt/latest mkl/2023.0.0 mpi/latest
    oclfpga/latest compiler/latest
NSIMUL is 4
===== AutoNEB Job Starting =====
 ===== Make Initial NEB Guess =====
---- Fix Atoms below 0.0 in direction y ----
---- Set initial magmom for [] to [] ----
---- Fix Atoms below 0.0 in direction y ----
---- Set initial magmom for [] to [] ----
---- Fix Atoms below 0.0 in direction y ----
---- Set initial magmom for [] to [] ----
---- Fix Atoms below 0.0 in direction y ----
---- Set initial magmom for [] to [] ----
--- Successfully make guessed image chain by idpp method ! ---
===== Running AutoNEB =====
Notice: AutoNEB method is set
You manually set n_simul = 4, n_max = 10
----- Running AutoNEB -----
----- improvedtangent method is being used -----
The NEB initially has 6 images  (including the end-points)
Start of evaluation of the initial images
Now starting iteration 1 on  [0, 1, 2, 3, 4, 5]
Finished initialisation phase.
****Now adding another image until n_max is reached (6/10)****
Adding image between 2 and 3. New image point is selected on the basis of the biggest spring length!
Now starting iteration 2 on  [1, 2, 3, 4, 5, 6]
****Now adding another image until n_max is reached (7/10)****
Adding image between 1 and 2. New image point is selected on the basis of the biggest spring length!
Now starting iteration 3 on  [1, 2, 3, 4, 5, 6]
****Now adding another image until n_max is reached (8/10)****
Adding image between 5 and 6. New image point is selected on the basis of the biggest spring length!
Now starting iteration 4 on  [2, 3, 4, 5, 6, 7]
****Now adding another image until n_max is reached (9/10)****
Adding image between 7 and 8. New image point is selected on the basis of the biggest spring length!
Now starting iteration 5 on  [4, 5, 6, 7, 8, 9]
n_max images has been reached
****Now doing the CI-NEB calculation****
Now starting iteration 6 on  [2, 3, 4, 5, 6, 7]
----- AutoNEB calculation finished -----
===== AutoNEB Process Done ! =====
===== Running Post-Processing =====
=== n_max set to 0, automatically detect the images of chain by NEBTools ===
Appears to be only one band in the images.
num: 0; Energy: -11866.824886 (eV)
num: 1; Energy: -11866.82722 (eV)
num: 2; Energy: -11866.801706 (eV)
num: 3; Energy: -11866.78262 (eV)
num: 4; Energy: -11866.587954 (eV)
num: 5; Energy: -11865.497 (eV)
num: 6; Energy: -11866.330196 (eV)
num: 7; Energy: -11866.396832 (eV)
num: 8; Energy: -11866.429487 (eV)
num: 9; Energy: -11866.435419 (eV)
Reaction Barrier and Energy Difference: (1.3278860000009587, 0.3894670000008773) (eV)
Appears to be only one band in the images.
Processing band          0 /          1
===== Done ! =====

      Step     Time          Energy          fmax
FIRE:    0 00:35:21   -11864.578532         3.642360
FIRE:    1 00:36:10   -11864.782954         2.870872
FIRE:    2 00:36:56   -11865.029191         2.407401
FIRE:    3 00:37:40   -11865.183575         3.420662
FIRE:    4 00:38:26   -11865.149576         4.842759
FIRE:    5 00:39:11   -11865.185455         3.786428
FIRE:    6 00:40:01   -11865.238631         2.154321
FIRE:    7 00:40:47   -11865.284643         1.979722
FIRE:    8 00:41:33   -11865.307675         1.889517
FIRE:    9 00:42:20   -11865.309583         1.742232
FIRE:   10 00:43:06   -11865.304529         1.941274
FIRE:   11 00:43:51   -11865.305786         1.824236
FIRE:   12 00:44:39   -11865.321340         1.566756
FIRE:   13 00:45:26   -11865.351452         1.232224
FIRE:   14 00:46:12   -11865.385412         1.163525
FIRE:   15 00:47:00   -11865.405953         1.863983
FIRE:   16 00:47:46   -11865.411933         1.196430
FIRE:   17 00:48:34   -11865.426725         1.069129
FIRE:   18 00:49:19   -11865.431487         0.961849
FIRE:   19 00:50:06   -11865.439472         0.760762
FIRE:   20 00:50:51   -11865.448191         0.517175
FIRE:   21 00:51:39   -11865.455199         0.459462
FIRE:   22 00:52:25   -11865.459289         0.472398
FIRE:   23 00:53:13   -11865.460873         0.431848
FIRE:   24 00:54:00   -11865.461725         0.539257
FIRE:   25 00:54:45   -11865.463809         0.619627
FIRE:   26 00:55:31   -11865.468630         0.598042
FIRE:   27 00:56:20   -11865.476511         0.468127
FIRE:   28 00:57:07   -11865.485968         0.365569
FIRE:   29 00:57:53   -11865.494549         0.412251
FIRE:   30 00:58:43   -11865.501988         0.385664
FIRE:   31 00:59:30   -11865.511484         0.422284
FIRE:   32 01:00:16   -11865.525009         0.366845
FIRE:   33 01:01:02   -11865.539757         0.231356
FIRE:   34 01:01:50   -11865.554350         0.309826
FIRE:   35 01:02:39   -11865.570950         0.294388
FIRE:   36 01:03:24   -11865.590207         0.275131
FIRE:   37 01:04:08   -11865.611798         0.247073
FIRE:   38 01:04:53   -11865.636943         0.299314
FIRE:   39 01:05:39   -11865.672444         0.295211
FIRE:   40 01:06:26   -11865.724616         0.398352
FIRE:   41 01:07:17   -11865.809261         0.672083
FIRE:   42 01:08:02   -11865.836765         0.677076
FIRE:   43 01:08:49   -11865.840971         0.681209
FIRE:   44 01:09:39   -11865.848278         0.688700
FIRE:   45 01:10:27   -11865.858156         0.697759
FIRE:   46 01:11:13   -11865.871159         0.704303
FIRE:   47 01:11:59   -11865.887878         0.703259
FIRE:   48 01:12:47   -11865.907860         0.688451
FIRE:   49 01:13:34   -11865.929898         0.651860
FIRE:   50 01:14:24   -11865.955356         0.574260
FIRE:   51 01:15:10   -11865.982482         0.436765
FIRE:   52 01:15:55   -11866.006105         0.314752
FIRE:   53 01:16:44   -11866.020034         0.413740
FIRE:   54 01:17:31   -11866.025234         0.697364
FIRE:   55 01:18:18   -11866.027275         0.653814
FIRE:   56 01:19:01   -11866.030891         0.572655
FIRE:   57 01:19:46   -11866.035293         0.464302
FIRE:   58 01:20:32   -11866.039733         0.341400
FIRE:   59 01:21:23   -11866.043688         0.219902
FIRE:   60 01:22:07   -11866.046996         0.200859
FIRE:   61 01:22:55   -11866.049775         0.236517
FIRE:   62 01:23:41   -11866.052505         0.255194
FIRE:   63 01:24:30   -11866.055451         0.245242
FIRE:   64 01:25:18   -11866.058968         0.220082
FIRE:   65 01:26:06   -11866.063303         0.193788

      Step     Time          Energy          fmax
FIRE:    0 01:57:54   -11866.084572         0.652487
FIRE:    1 01:58:40   -11866.076705         0.676456
FIRE:    2 01:59:28   -11866.059007         0.715741
FIRE:    3 02:00:13   -11866.028561         0.745702
FIRE:    4 02:00:59   -11865.982395         0.745921
FIRE:    5 02:01:46   -11865.918309         0.713379
FIRE:    6 02:02:36   -11865.835733         0.659343
FIRE:    7 02:03:21   -11865.736940         0.602895
FIRE:    8 02:04:09   -11865.616163         0.570285
FIRE:    9 02:04:56   -11865.486295         0.576872
FIRE:   10 02:05:41   -11865.366292         0.553076
FIRE:   11 02:06:26   -11865.290706         0.491710
FIRE:   12 02:07:10   -11865.268212         0.382165
FIRE:   13 02:07:56   -11865.295327         0.349660
FIRE:   14 02:08:42   -11865.381779         0.505350
FIRE:   15 02:09:29   -11865.389766         0.478347
FIRE:   16 02:10:14   -11865.402261         0.431838
FIRE:   17 02:11:01   -11865.414176         0.381902
FIRE:   18 02:11:46   -11865.421392         0.382459
FIRE:   19 02:12:33   -11865.422792         0.386433
FIRE:   20 02:13:18   -11865.419163         0.378440
FIRE:   21 02:14:05   -11865.411475         0.372510
FIRE:   22 02:14:52   -11865.399699         0.363692
FIRE:   23 02:15:36   -11865.382875         0.311966
FIRE:   24 02:16:22   -11865.360650         0.348324
FIRE:   25 02:17:08   -11865.335546         0.494877
FIRE:   26 02:17:55   -11865.316262         0.735960
FIRE:   27 02:18:45   -11865.307194         1.100647
FIRE:   28 02:19:30   -11865.309266         0.196186
FIRE:   29 02:20:19   -11865.313339         0.909086
FIRE:   30 02:21:08   -11865.313751         0.707882
FIRE:   31 02:21:54   -11865.314565         0.350729
FIRE:   32 02:22:42   -11865.315812         0.178845
FIRE:   33 02:23:30   -11865.315906         0.178576
FIRE:   34 02:24:16   -11865.316091         0.178013
FIRE:   35 02:25:07   -11865.316373         0.177024
FIRE:   36 02:25:57   -11865.316741         0.175838
FIRE:   37 02:26:45   -11865.317208         0.174313
FIRE:   38 02:27:32   -11865.317762         0.172504
FIRE:   39 02:28:19   -11865.318410         0.170232
FIRE:   40 02:29:08   -11865.319231         0.167451
FIRE:   41 02:29:54   -11865.320258         0.163917
FIRE:   42 02:30:40   -11865.321550         0.165854
FIRE:   43 02:31:28   -11865.323145         0.171715
FIRE:   44 02:32:18   -11865.325136         0.179361
FIRE:   45 02:33:07   -11865.327600         0.188798
FIRE:   46 02:33:54   -11865.330675         0.200416
FIRE:   47 02:34:43   -11865.334504         0.213338
FIRE:   48 02:35:32   -11865.339264         0.226164
FIRE:   49 02:36:21   -11865.345094         0.234130
FIRE:   50 02:37:09   -11865.352065         0.232614
FIRE:   51 02:37:54   -11865.360183         0.216983
FIRE:   52 02:38:43   -11865.369363         0.185817
FIRE:   53 02:39:29   -11865.379371         0.148843
FIRE:   54 02:40:14   -11865.389614         0.138991
FIRE:   55 02:40:59   -11865.399380         0.137180
FIRE:   56 02:41:45   -11865.408222         0.165171
FIRE:   57 02:42:32   -11865.416086         0.203523
FIRE:   58 02:43:19   -11865.423394         0.231605
FIRE:   59 02:44:06   -11865.430772         0.247774
FIRE:   60 02:44:53   -11865.438325         0.254307
FIRE:   61 02:45:41   -11865.445080         0.267498
FIRE:   62 02:46:29   -11865.448719         0.881362
FIRE:   63 02:47:15   -11865.451780         0.231961
FIRE:   64 02:48:02   -11865.452007         0.231221
FIRE:   65 02:48:48   -11865.452438         0.229862
FIRE:   66 02:49:34   -11865.453071         0.227879
FIRE:   67 02:50:21   -11865.453876         0.225257
FIRE:   68 02:51:08   -11865.454816         0.222156
FIRE:   69 02:52:00   -11865.455865         0.218612
FIRE:   70 02:52:44   -11865.457001         0.214634
FIRE:   71 02:53:30   -11865.458341         0.209595
FIRE:   72 02:54:15   -11865.459892         0.203370
FIRE:   73 02:55:01   -11865.461651         0.195673
FIRE:   74 02:55:49   -11865.463637         0.186376
FIRE:   75 02:56:36   -11865.465861         0.175254
FIRE:   76 02:57:22   -11865.468343         0.162491
FIRE:   77 02:58:09   -11865.471156         0.147879
FIRE:   78 02:58:58   -11865.474326         0.131731
FIRE:   79 02:59:44   -11865.477956         0.114483
FIRE:   80 03:00:31   -11865.482135         0.096221
FIRE:   81 03:01:17   -11865.486787         0.084464
FIRE:   82 03:02:05   -11865.491684         0.134461
FIRE:   83 03:02:50   -11865.496380         0.273588
FIRE:   84 03:03:36   -11865.496522         0.051535
FIRE:   85 03:04:21   -11865.496806         0.203176
FIRE:   86 03:05:08   -11865.496838         0.161525
FIRE:   87 03:05:53   -11865.496901         0.086825
FIRE:   88 03:06:42   -11865.497000         0.040039

# JamesMisaka in 2023-11-30
# Vibrational analysis from finite displacement by using abacus
# part of ATST-Tools scripts

import os
import numpy as np
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo
from ase.io import read, write
from ase.calculators.abacus import Abacus, AbacusProfile
from ase.parallel import world, parprint

stru = "STRU"
atoms = read(stru)
# indices setting for which atoms to be displaced
#vib_indices = [atom.index for atom in atoms if atom.symbol == 'H']
vib_indices = [0, 1, 37]
T = 523.15 # K

vib_name = 'vib'
delta = 0.01
nfree = 2

abacus = "abacus"
mpi = 16
omp = 4
lib_dir = "/lustre/home/2201110432/example/abacus"
pseudo_dir = f"{lib_dir}/PP"
basis_dir = f"{lib_dir}/ORB"
pp = {
      'H': 'H_ONCV_PBE-1.0.upf',
      'C': 'C_ONCV_PBE-1.0.upf',
      'O': 'O_ONCV_PBE-1.0.upf',
      'Fe': 'Fe_ONCV_PBE-1.0.upf',
      }
basis = {
         'H': 'H_gga_6au_100Ry_2s1p.orb',
         'C': 'C_gga_7au_100Ry_2s2p1d.orb',
         'O': 'O_gga_7au_100Ry_2s2p1d.orb',
         'Fe': 'Fe_gga_8au_100Ry_4s2p2d1f.orb'
         ,}
kpts = [3, 1, 2]
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'none',
    'smearing_method': 'mp',
    'smearing_sigma': 0.002,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'scf_thr': 1e-7,
    'scf_nmax': 300,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'init_chg': 'atomic',
    'init_wfc': 'file',
    'cal_force': 1,
    'cal_stress': 1,
    'out_stru': 1,
    'out_chg': 1,
    'out_mul': 0,
    'out_bandgap': 0,
    'out_wfc_lcao': 1,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
    'efield_pos_max': 0.7
}


def set_calculator(abacus, parameters, mpi=1, omp=1) -> Abacus:
    """Set Abacus calculators"""
    os.environ['OMP_NUM_THREADS'] = f'{omp}'
    profile = AbacusProfile(
        argv=['mpirun', '-np', f'{mpi}', abacus])
    out_directory = f"SCF-rank{world.rank}"
    calc = Abacus(profile=profile, directory=out_directory,
                **parameters)
    return calc

if __name__ == "__main__":
    print("==> Starting Vibrational Analysis <==")

    atoms.calc = set_calculator(abacus, parameters, mpi=mpi, omp=omp)

    vib = Vibrations(atoms, indices=vib_indices, 
                    name=vib_name, delta=delta, nfree=nfree)

    print("==> Running Vibrational Analysis <==")
    vib.run()
    # post-processing
    print("==> Done !!! <==")
    print(f"==> All force cache will be in {vib_name} directory <==")
    print("==> Vibrational Analysis Summary <==")
    vib.summary()
    print("==> Writing All Mode Trajectory <==")
    vib.write_mode()
    # thermochemistry
    print("==> Doing Harmonmic Thermodynamic Analysis <==")
    real_vib_energies = np.array([energy for energy in vib.get_energies() 
                        if energy.imag == 0 and energy.real > 0], dtype=float)
    thermo = HarmonicThermo(real_vib_energies)
    entropy = thermo.get_entropy(T)
    free_energy = thermo.get_helmholtz_energy(T)
    print(f"==> Entropy: {entropy:.6e} eV/K <==")
    print(f"==> Free Energy: {free_energy:.6f} eV <==")
    

Loading mkl version 2023.0.0
Loading tbb version 2021.8.0
Loading compiler-rt version 2023.0.0
Loading mpi version 2021.8.0
Loading compiler version 2023.0.0
Loading oclfpga version 2023.0.0
  Load "debugger" to debug DPC++ applications with the gdb-oneapi debugger.
  Load "dpl" for additional DPC++ APIs: https://github.com/oneapi-src/oneDPL

Loading abacus/3.4.3-icx-dev
  Loading requirement: tbb/latest compiler-rt/latest mkl/latest mpi/latest
    oclfpga/latest compiler/latest
Vibrational Calculation Start !
==> Running Vibrational Analysis <==
==> Get Energy and Frequency <==
==> Get ZPE <==
==> Writing Trajectory <==
==> Summary <==
---------------------
  #    meV     cm^-1
---------------------
  0   87.4i    705.0i
  1    2.4      19.3
  2    3.5      28.3
  3    6.1      48.9
  4   14.5     117.2
  5   17.4     140.7
  6   20.9     168.5
  7   30.4     244.9
  8   38.7     312.0
  9   51.0     411.6
 10   53.1     428.5
 11   57.8     466.1
 12   65.3     526.4
 13   79.5     641.3
 14   96.8     781.1
 15   98.2     791.9
 16  102.3     825.1
 17  107.3     865.6
 18  108.3     873.4
 19  111.7     900.6
 20  120.9     975.1
 21  123.8     998.4
 22  129.7    1045.8
 23  131.5    1060.4
 24  132.7    1070.3
 25  137.3    1107.4
 26  144.5    1165.2
 27  151.3    1220.2
 28  153.9    1241.2
 29  155.6    1255.2
 30  158.1    1275.4
 31  159.4    1285.7
 32  163.2    1316.1
 33  163.7    1320.7
 34  165.8    1337.5
 35  166.6    1343.8
 36  167.1    1348.0
 37  178.3    1438.0
 38  178.9    1442.9
 39  179.2    1445.6
 40  179.8    1450.5
 41  181.3    1462.1
 42  263.0    2121.5
 43  362.0    2920.0
 44  366.9    2959.6
 45  367.8    2966.6
 46  368.1    2968.8
 47  369.0    2976.1
 48  371.0    2992.1
 49  373.5    3012.7
 50  374.4    3019.4
 51  374.7    3022.4
 52  375.8    3031.3
 53  377.5    3044.9
---------------------
Zero-point energy: 4.416 eV
===== Done at Wed Nov 29 23:13:33 CST 2023! =====