Merge pull request #5 from QuantumMisaka/neb-dimer

Neb-dimer workflow

Author: QuantumMisaka

README.md

@@ -1,8 +1,11 @@
-# ase-neb-abacus
-NEB or CI-NEB workflow by ASE and ABACUS using ASE-ABACUS interface,
-Which use enhanced NEB method in ASE, like Dynamic NEB `DyNEB`, improved-tangent method `IT-NEB`, and also `AutoNEB` in ASE.  
+# ATST-Tools
+Advanced ASE Transition State Tools for ABACUS, including:
+- NEB and improvement, like Dynamic NEB
+- CI-NEB, IT-NEB and others
+- AutoNEB
+- Dimer
 
-Version v1.0.2
+Version v1.1
 
 ## Dependencies:
 - [ASE](https://wiki.fysik.dtu.dk/ase/about.html)
@@ -14,12 +17,13 @@ Notice: GPAW and ABACUS should be dependent on same MPI and libraries environmen
 For instance, if your ABACUS is installed by Intel-OneAPI toolchain, your GPAW should NOT be dependent on gcc-toolchain like OpenMPI and OpenBLAS.
 
 
-## Usage
-The workflow is based on 3 main python scripts and 1 workflow submit script. Namely:
+## NEB workflow
+### Usage
+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, 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}` for each rank which do calculation of each image. 
+- `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.
 
@@ -43,13 +47,25 @@ Also, user can run each step in one script `neb_submit.sh` by `bash neb_submit.s
 
 > Notice: Before you start neb calculation process, make sure that you have check the nodes and cpus setting in `neb_submit.sh` and `neb_run.py` to make sure that you can reach the highest performance !!! 
 
-## Method
+### 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, and `ase.mep.autoneb.AutoNEB` method should be independently 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. But in `AutoNEB`, `eb` method is used for default.
 - 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: 
 - 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.
 
+## Dimer workflow
+### 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
+- - `displacement_vector.npy` for displacement vector of Dimer calculation  
+- `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.
+- `neb-dimer_srun.sh` is a try to run NEB + Dimer calculation in one server scripts. The Default setting is for `slurm` job submission system.
+
+### Method
+(Waiting for update)
 
 ## Examples
 - Li-diffu-Si: Li diffusion in Si, an example for running ASE-NEB-ABACUS based on existing ABACUS input files of initial and final state, using ABACUS as SCF calculator and ASE as optimizer and NEB calculator.  Also, an dflow example is proposed.
@@ -58,9 +74,9 @@ Also, user can run each step in one script `neb_submit.sh` by `bash neb_submit.s
 - CO-Pt111 : CO dissociation on Pt(111) surface, an example for running ASE-NEB-ABACUS based on existing ABACUS input files of initial and final state, use ABACUS as SCF calculator, use ABACUS as optimizer for optimization of initial state and final state, use ASE as NEB calculator
 - Cy-Pt_graphene: Cyclohexane dehydrogenation on Pt-doped graphene surface, an example for running ASE-NEB-ABACUS based on existing ABACUS input files of initial and final state, use ABACUS as optimizer for optimization of initial state and final state, use ASE as NEB calculator
 
+AutoNEB example ix on update, Dimer example is preparing
 
 ## Next Examples
-- `AutoNEB` test 
 - CO-Fe100: CO dissociation on Fe(100) surface, in this example we will focus on setting `magmom` during NEB calculation, which is important for spin-polarized magnetic system.
 - TS of MTM process by H2O2 in Cu/Ag-ZSM5 system, including:
 - - proton transfer process
@@ -83,7 +99,8 @@ Also, user can run each step in one script `neb_submit.sh` by `bash neb_submit.s
 - [x] More test in surface reaction system
 - [x] Parallel computing during images relaxation by `gpaw python`
 - [x] `AutoNEB` implementation
+- [x] Connected to Dimer method
 - [ ] More test in magnetic surface reaction system
-- [ ] Connected to Dimer method
+- [ ] put calculation setting in an independent file (decoupling run*.py)
 
 

dimer/dimer_run.py

@@ -0,0 +1,155 @@
+from ase import Atom, Atoms
+from ase.io import Trajectory, read, write
+from ase.mep import DimerControl, MinModeAtoms, MinModeTranslate
+from ase.optimize import BFGS, FIRE
+from ase.calculators.abacus import Abacus, AbacusProfile
+import os, sys
+
+import numpy as np
+
+fmax = 0.05
+dimer_input_file = 'dimer_init.traj'
+init_eigenmode_method = 'displacement'
+displacement_input = 'displacement_vector.npy'
+
+# setting for calculator
+mpi = 16
+omp = 4
+abacus = "abacus"
+#lib_dir = "/lustre/home/2201110432/example/abacus"
+lib_dir = "/home/james/example"
+pseudo_dir = f"{lib_dir}/PP"
+basis_dir = f"{lib_dir}/ORB"
+pp = {
+        'H':'H_ONCV_PBE-1.0.upf',
+        'Au':'Au_ONCV_PBE-1.0.upf',
+      }
+basis = {
+        'H':'H_gga_6au_100Ry_2s1p.orb',
+        'Au':'Au_gga_7au_100Ry_4s2p2d1f.orb',
+        }
+kpts = [3, 1, 3]
+parameters = {
+    'calculation': 'scf',
+    'nspin': 2,
+    'xc': 'pbe',
+    'ecutwfc': 100,
+    '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,
+    'out_stru': 1,
+    'out_chg': 0,
+    'out_mul': 0,
+    'out_bandgap': 0,
+    'efield_flag': 1,
+    'dip_cor_flag': 1,
+    'efield_dir': 1,
+    'efield_pos_max': 0.7
+}
+
+
+class AbacusDimer:
+    """Customize Dimer calculation by using ABACUS"""
+    
+    def __init__(self, init_Atoms, parameters, abacus='abacus',
+                 mpi=1, omp=1, directory='DIMER', 
+                 traj_file='dimer.traj',
+                 init_eigenmode_method='displacement',
+                 displacement_vector: np.ndarray = None,):
+        """Initialize Dimer method by using ASE-ABACUS
+
+        init_Atoms (Atoms object): starting image, can be from every way including NEB result
+        parameters (dict): settings of abacus input parameters
+        abacus (str): Abacus executable file. Default: 'abacus'
+        directory (str): calculator directory name, for parallel calculation {directory}-rank{i} will be the directory name
+        mpi (int): number of MPI for abacus calculator
+        omp (int): number of OpenMP for abacus calculator
+        traj_file (str): trajectory file name for dimer calculation, when running dimer calculation, trajetory will be written to this file, default is 'dimer.traj'
+        init_eigenmode_method (str): dimer initial eigenmode method. Choose from 'displacement' and 'gauss'.
+        displacement_vector (np.ndarray): displacement vector for dimer initial eigenmode. Only used when init_eigenmode_method is 'displacement'
+        """
+        self.init_Atoms = init_Atoms
+        self.parameters = parameters
+        self.abacus = abacus
+        self.mpi = mpi
+        self.omp = omp
+        self.directory = directory
+        self.traj_file = traj_file
+        self.init_eigenmode_method = init_eigenmode_method
+        self.displacement_vector = displacement_vector
+        
+    def set_calculator(self):
+        """Set Abacus calculators"""
+        os.environ['OMP_NUM_THREADS'] = f'{self.omp}'
+        profile = AbacusProfile(
+            argv=['mpirun', '-np', f'{self.mpi}', self.abacus])
+        out_directory = self.directory
+        calc = Abacus(profile=profile, directory=out_directory,
+                        **self.parameters)
+        return calc
+    
+    def set_d_mask_by_displacement(self):
+        """set mask by displacement"""
+        print("=== Set mask by displacement vector where displacement is [0,0,0] ===")
+        d_mask = self.displacement_vector != np.zeros(3)
+        d_mask = d_mask[:,0].tolist()
+        return d_mask
+    
+    def set_d_mask_by_constraint(self):
+        """set mask by constraint of Atoms"""
+        print("=== Set mask by constraint read from init Atoms ===")
+        dimer_init = self.init_Atoms
+        const_list = dimer_init._get_constraints()[0]
+        const_list = const_list.todict()['kwargs']['indices']
+        d_mask = [True] * len(dimer_init)
+        for ind in const_list:
+            d_mask[ind] = False
+        return d_mask
+        
+    def run(self, fmax=0.05):
+        """run dimer calculation workflow"""
+        dimer_init = self.init_Atoms
+        dimer_init.calc = self.set_calculator()
+        dimer_traj = Trajectory(self.traj_file, 'w', dimer_init)
+        if self.init_eigenmode_method == "displacement":
+            # d_mask = self.set_d_mask_by_displacement() # not the best
+            d_mask = self.set_d_mask_by_constraint()
+            d_control = DimerControl(initial_eigenmode_method=self.init_eigenmode_method, 
+                                    displacement_method="vector", 
+                                    mask=d_mask)
+            d_atoms = MinModeAtoms(dimer_init, d_control)
+            d_atoms.displace(displacement_vector=self.displacement_vector)
+        elif self.init_eigenmode_method == "gauss":
+            # leave a way for random displacement
+            d_mask = self.set_d_mask_by_constraint()
+            d_control = DimerControl(initial_eigenmode_method=self.init_eigenmode_method, 
+                                    mask=d_mask)
+            d_atoms = MinModeAtoms(dimer_init, d_control)
+        else:
+            raise ValueError("init_eigenmode_method must be displacement or gauss")
+        dimer_relax = MinModeTranslate(d_atoms, trajectory=dimer_traj)
+        dimer_relax.run(fmax=fmax)
+        
+if __name__ == "__main__":
+# running process
+# read initial guessed neb chain
+    dimer_init = read(dimer_input_file)
+    displacement_vector = np.load(displacement_input)
+    dimer = AbacusDimer(dimer_init, parameters, abacus=abacus,
+                        mpi=mpi, omp=omp, 
+                        init_eigenmode_method=init_eigenmode_method,
+                        displacement_vector=displacement_vector)
+    dimer.run(fmax=fmax)

dimer/dimer_submit.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=16
+#SBATCH --cpus-per-task=4
+#SBATCH -J DIMER-ABACUS
+#SBATCH -o run_dimer.out
+#SBATCH -e run_dimer.err
+#SBATCH -p C064M0256G
+
+# JamesMisaka in 2023-11-14
+# workflow of ase-abacus-dimer method
+# for one calculator, ntasks-per-node for mpi, cpus-per-task for openmp
+
+# in developer's PKU-WM2 server
+source /lustre/home/2201110432/apps/miniconda3/etc/profile.d/conda.sh
+conda activate gpaw-intel
+module load abacus/3.4.2-icx
+
+# if one just done neb calculation
+python neb2dimer.py neb.traj # or neb_latest.traj
+
+# Job state 
+echo $SLURM_JOB_ID > JobRun.state
+echo "Start at $(date)" >> JobRun.state
+
+# just do dimer by run !
+# do check your ntasks(mpi) and cpu(omp) setting is right
+python dimer_run.py
+
+echo "===== Dimer Calculation Done at $(date)! ====="
+
+# Job State
+echo "End at $(date)" >> JobRun.state

dimer/neb-dimer_srun.sh

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+# full workflow controling NEB and Dimer calculation
+
+NEB_JOBID=$(sbatch neb_submit.sh) | awk '{print $4}' # get the jobid of neb and submit neb job
+echo "NEB_JOBID is ${NEB_JOBID}"
+DIMER_JOBID=$(sbatch --dependency=afterok:${NEB_JOBID} dimer_submit.sh) | awk '{print $4}' # submit dimer job after neb job done
+echo "Dimer Job Submitted and will be run after NEB job !"
+echo "DIMER_JOBID is ${DIMER_JOBID}"

dimer/neb2dimer.py

@@ -0,0 +1,73 @@
+import numpy as np
+import sys
+from ase.io import Trajectory, read, write
+from ase.mep.neb import NEBTools
+
+def neb2dimer(neb_traj:list, n_max:int=0, out_traj='dimer_init.traj', out_vec='displacement_vector.npy', 
+                step_before_TS:int=1, step_after_TS:int=1):
+    '''
+    Transform neb chain to dimer init
+    '''
+    # you should use the latest neb chain or other single neb chain
+    # but a full neb chain is permitted, in which the latest neb chain will be used
+    if n_max == 0:
+        print("=== n_max set to 0, get n_images by using NEBTools ===")
+        n_images = NEBTools(neb_traj)._guess_nimages()
+    elif (n_max >= 0) and (type(n_max) == int):
+        n_images = n_max + 2
+    else:
+        raise ValueError("n_max must be a non-negative integer")
+    if len(neb_traj) < n_images:
+        raise ValueError(f"n_images ={n_images} is larger than the length of neb_traj ={len(neb_traj)}")
+    elif len(neb_traj) > n_images:
+        print(f"=== n_images ={n_images} is smaller than the length of neb_traj ={len(neb_traj)} ! Only the last {n_images} images are used ===")
+    # used neb_chain is the final traj
+    neb_chain = neb_traj[-n_images:]
+    # get TS information from NEB chain
+    barrier = NEBTools(neb_chain).get_barrier()[0]
+    fmax = NEBTools(neb_chain).get_fmax()
+    raw_barrier = max([image.get_potential_energy() for image in neb_chain])
+    TS_info = [(ind, image) 
+            for ind, image in enumerate(neb_chain) 
+            if image.get_potential_energy() == raw_barrier][0]
+    print(f"=== Locate TS in {TS_info[0]} of 0-{n_images-1} images  ===")
+    print(f"=== TS Barrier: {barrier:.4f} (eV) ===")
+    print(f"=== TS Fmax: {fmax:.4f} (eV/A) ===")
+    print(f"=== TS images is output as {out_traj} for dimer init ===")
+    # output TS of neb for dimer init
+    write(out_traj, TS_info[1], format='traj')
+    # output displancement vector by using the nearest two images of TS
+    ind_before_TS = TS_info[0] - step_before_TS
+    ind_after_TS = TS_info[0] + step_after_TS
+    img_before = neb_chain[ind_before_TS]
+    img_after = neb_chain[ind_after_TS]
+    print(f"=== Displacement vector is generated by position minus from {ind_after_TS} image to {ind_before_TS} image ===")
+    displacement_vector = (img_after.positions - img_before.positions)
+    print(f"=== Displacement vector is output as {out_vec} for dimer init ===")
+    np.save(out_vec,displacement_vector)
+    return
+    
+if __name__ == "__main__":
+    msg = '''
+neb2dimer.py is to make dimer inputfile by using neb result
+Usage: 
+    python neb2dimer.py [traj_file] ([n_max])
+        Notice that n_max = n_images - 2
+These sctipt will output two files:
+    dimer_init.traj: TS of neb chain for dimer init-structure
+    displacement_vector.npy: displacement vector of dimer init
+    '''
+    if len(sys.argv) < 2:
+        print(msg)
+    elif len(sys.argv) == 2:
+        if sys.argv == "--help" or sys.argv == "-h":
+            print(msg)
+        else:
+            traj_file = sys.argv[1]
+            neb_traj = Trajectory(traj_file)
+            neb2dimer(neb_traj)
+    else:
+        traj_file = sys.argv[1]
+        n_max = int(sys.argv[2])
+        neb_traj = Trajectory(traj_file)
+        neb2dimer(neb_traj, n_max)