ASE-ABACUS | 第二章：NEB过渡态计算与ATST-Tools

# ATST-Tools
Advanced ASE Transition State Tools for ABACUS and Deep-Potential, including:
- NEB, including CI-NEB, IT-NEB and others.
- Serial Dynamic NEB (DyNEB) calculation.
- AutoNEB: an automatic NEB workflow.
- Single-End TS search: Sella, Dimer.
- IRC analysis by Sella.
- Double-to-single (D2S) TS workflow: neb2dimer, neb2sella.
- Vibration analysis and ideal gas thermochemistry analysis.

Version 1.5.0

Copyright @ QuantumMisaka from TMC-PKU & AISI

## Update log from 1.4.0 to 1.5.0
- Add internal reaction coordination (IRC) calculation by using Sella packages.
- Add *neb2dimer_dp.py*, *neb2sella_dp.py*, *relax_dp.py* for deepmd usage, including tf and torch version, and all D2S method are using pymatgen IDPP.
- Change default relaxation optimizer from BFGS to QuasiNewton (BFGSLineSearch).
- change the way to use `AbacusProfile` to march the newest version of ase-abacus interface, which use different command format.
- Since constraints read-in problem in ASE-ABACUS interface is fixed [#28](https://gitlab.com/1041176461/ase-abacus/-/issues/28), the constraints is the two scripts above is modified, make Dimer and Sella useful.
- Single-End TS search scripts update using Dimer method and Sella packages.
- Default *neb_make.py* change to pymatgen version below.
- *neb_make_pymatgen.py* scripts by @MoseyQAQ, which use IDPP of pymatgen to do NEB initial guess, avoiding IDPP problem from ASE including edge-crossing problem in NEB guess generation.
- Support DeepPotential and DPA-2 Potential usage scripts in `ase-dp` directory.
- New function: `neb2vib` method to use NEB chain information to do partial vibrational analysis.
- Thermochemistry analysis for ideal gas system.
- Examples Update

## Dependencies:
- [ASE](https://wiki.fysik.dtu.dk/ase/about.html), but you should install ASE by [ASE-ABACUS interface](https://gitlab.com/1041176461/ase-abacus) as a separate ASE package for ABACUS usage.
- [ABACUS](https://github.com/deepmodeling/abacus-develop), one can install ABACUS by ABACUS toolchain [ABACUS](https://github.com/deepmodeling/abacus-develop/tree/develop/toolchain), or refer to [ABACUS-docs](https://abacus.deepmodeling.com/en/latest/)
- [pymatgen](https://pymatgen.org/) and [pymatgen-analysis-diffusion](https://github.com/materialsvirtuallab/pymatgen-analysis-diffusion) in the usage of new `neb_make.py` script and D2S TS method, which can be installed by `pip install pymatgen pymatgen-analysis-diffusion`
- [Sella](https://github.com/zadorlab/sella) if one wants to use Sella method for Single-End TS search. which can be installed by `pip install sella`
- [GPAW](https://wiki.fysik.dtu.dk/gpaw/install.html) if one wants to run NEB images relaxation in parallel. The installation of GPAW need some efferts, and one can refer to [GPAW installation](https://wiki.fysik.dtu.dk/gpaw/install.html) for more details.
- [deepmd-kit](https://github.com/deepmodeling/deepmd-kit) or [DPA-2](https://zenodo.org/records/10428497) if one wants to use Deep-Potential or DPA-2 potential, one can refer to [deepmd-docs](https://docs.deepmodeling.com/projects/deepmd/en/master/)


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.

## Tutorials
- One can get ASE-ABACUS usage from this [Bohrium Notebook 1](https://bohrium.dp.tech/notebooks/6516485694)
- One can get NEB tutorial from this [Bohrium Notebook 2](https://nb.bohrium.dp.tech/detail/39369325971)
- One can get Single-End TS search tutorial from this [Bohrium Notebook 3](https://bohrium.dp.tech/notebooks/29581597682)
- D2S TS exploration tutorial is coming soon.

## Workflow

![ATST-NEB-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
└── neb2vib.py
```
Before use running scripts, you should add these libraries into your PYTHONPATH:
```bash
export PYTHONPATH=/path/to/source:$PYTHONPATH
```

There are also other workflow in ATST-Tools
- AutoNEB TS exploration
- Double-to-Single (D2S) TS exploration

## Developing
- [ ] 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.
- [ ] Move workflow parts of Deep-Potential and DPA-2 potential to `source` directory
- [ ] Parallel NEB calculation for D2S method and DP usage
- [ ] [dflow](https://github.com/dptech-corp/dflow) version of ATST-Tools workflow
- [ ] Bond soft scanning and automatic reaction locating (proposed by [CLAM workflow](https://github.com/lalaheihaihei/catalyticLAM))
- [x] Test for Sella usage to get best performance
- [x] Checkout the problem in Dimer-ABACUS calculation
- [x] Inplement `neb2vib` method in  `neb2dimer` workflow for DPA-2 (and ABACUS)
- [x] Other TS method usage, like [Sella](https://github.com/zadorlab/sella)


## 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 -i [INIT/result] [FINAL/result] -n [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 [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.

## Single-End TS search

Contrary to NEB, the single-end TS search is based on the exploration of saddle point by using the gradient and Hessian (or only Hessian eigenmode) information of points in PES, which have good efficiency if one have an approximate TS information. ATST-Tools support Dimer method and Sella method.

### Sella workflow
The Sella workflow is based on 1 main python scripts and 1 submit script, namely:
- `sella_run.py`  is the key running script of Dimer calculation, which will run Sella calculation based on supplying `STRU` files for initial state of Sella calculation. This script will generate `run_sella.traj` for Sella calculation trajectory. Users should edit this file to set parameters for Sella calculation, and run Sella calculation by `python sella_run.py`. When running, any Sella images calculation will be doing in `ABACUS` directory.
- `sella_submit.sh` will do Sella workflow in one scripts. The Default setting is for `slurm` job submission system.

Also, Sella package can be used for IRC calculation, one can refer to `sella_IRC.py` for more details.

For more detail on Sella principle and usage, one can refer to [Sella](https://github.com/zadorlab/sella) and [Sella-wiki](https://github.com/zadorlab/sella/wiki)


### Dimer workflow

The Dimer workflow is based on 2 main python scripts and 1 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. 
- - `STRU` files, representing the initial state of Dimer calculation, stored as `STRU` files for ABACUS calculation.
- - `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 `run_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 `ABACUS` directory.
- `dimer_submit.sh` will do Dimer workflow in one scripts. The Default setting is for `slurm` job submission system.


## Double-to-single (D2S) TS exploration

ATST-Tools also offer a double-to-single (D2S) TS exploration method, which is named as `neb2dimer_abacus.py` and `neb2sella_abacus.py` in `dimer` and `sella` directories respectively, which run a rough NEB first and use the maximum information for initial guess of running single-ended method like dimer or Sella, obtaining a high-efficiency TS exploration. 

In D2S workflow, Dynamic NEB acceleration is used for NEB calculation for serial usage for faster rough NEB calculation.

This workflow is only for serial NEB calculation, but is efficient enough for TS exploration. The parallel NEB acceleration is in considering.

## 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. 
> Notes: QuasiNewton method in ASE is BFGSLineSearch, which is default usage in ATST-Tools and can have robust structural relaxation ability.

## ase-dp
There are also scripts for Deep-Potential and DPA-2 potential usage in `ase-dp` directory for TS exploration and structural relaxation, some newest examples can be used by `neb2dimer_dp.py`, `neb2sella_dp.py` and `relax_dp.py` for Deep-Potential and DPA-2 potential usage. 

Parallel NEB version for ase-dp is also in developing.

## Examples
- Li-diffu-Si: Li diffusion in Si, very easy example for serial and parallel NEB calculation about diffusion system. Note that ATS-Tools will tune the parameters towards TS exploration, for those who want to do diffusion calculation, one should check the effectiveness of the parameters by yourself.
- H2-Au111: H2 dissociation on Au(111) surface. which will have NEB, serial DyNEB, AutoNEB, Dimer, Sella and D2S example. The barrier is around 1.1 eV consistent with existing paper and calculation result. IRC results shows the reaction path is correct.
- CO-Pt111 : CO dissociation on Pt(111) surface. which have high barrier as 1.5 eV and including diffusion part. AutoNEB will always fail due to the low starting image. By performing Dimer / Sella calculation, or use the D2S method, the TS can be explored in an efficient and accurate way. By using Sella package, one can also get the IRC of C-O dissociaition reaction process.
- 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. Dimer calculation also get wrong result in this case, while Sella perform good result.

More examples is welcomed from users. 

## 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 trajectory file, one should get it by `get_force(apply_constraint=False)`.
- in Dimer calculation, the Energy, Forces and Stress property will be stored in trajetory file after specified which properties need to be stored. (Sella calculation should be likely) (The most easiest lost information is the stress information)
- 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.

# 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}/PP"
basis_dir = f"{lib_dir}/ORB"
# default pp and basis is supported by ase-abacus interface, need to check usage
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,
}


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 = "/lustre/home/2201110432/example/abacus"
pseudo_dir = f"{lib_dir}/PP"
basis_dir = f"{lib_dir}/ORB"
# default pp and basis is supported by ase-abacus interface, need to check usage
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': 'd3_bj',
    'smearing_method': 'mp',
    'smearing_sigma': 0.002,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'mixing_beta': 0.4,
    'mixing_gg0': 1.0,
    'mixing_ndim': 8,
    'scf_thr': 1e-7,
    'scf_nmax': 300,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'cal_force': 1,
    'cal_stress': 0,
    'init_wfc': 'atomic',
    'init_chg': 'atomic',
    'out_stru': 1,
    'out_chg': -1,
    'out_mul': 1,
    'out_wfc_lcao': 0,
    'out_bandgap': 1,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
}



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
from ase.mep.neb import NEBTools
from neb2vib import neb2vib

# reading usage can be changed by user
# usage by neb
neb_traj = read('neb_latest.traj', index=':')
atoms, vib_indices = neb2vib(neb_traj)

# traditional
# stru = "STRU"
# atoms = read(stru)
# vib_indices = [0, 1, 37]

# indices setting for which atoms to be displaced
# vib_indices = [atom.index for atom in atoms if atom.symbol == 'H']

T = 523.15 # K

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"
# default pp and basis is supported by ase-abacus interface
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',
    'mixing_ndim': 8,
    'scf_thr': 1e-7,
    'scf_nmax': 300,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'init_chg': 'file',
    'init_wfc': 'atomic',
    'cal_force': 1,
    'cal_stress': 1,
    'out_stru': 1,
    'out_chg': 1,
    'out_mul': 0,
    'out_bandgap': 0,
    'out_wfc_lcao': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
}

# developer only
vib_name = 'vib'
delta = 0.01
nfree = 2

def set_calculator(abacus, parameters, mpi=1, omp=1) -> Abacus:
    """Set Abacus calculators"""
    os.environ['OMP_NUM_THREADS'] = f'{omp}'
    profile = AbacusProfile(f"mpirun -np {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 <==")
    vib_energies = vib.get_energies()
    #real_vib_energies = np.array([energy for energy in vib.get_energies() if energy.imag == 0 and energy.real > 0], dtype=float)
    thermo = HarmonicThermo(vib_energies, ignore_imag_modes=True,)
    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! =====