1+ '''this example shows how to use the IDPP method interpolate the
2+ nudged elastic band (NEB) calculations such that the transition
3+ state can be founded more quickly.
4+
5+ In this example, we also use the reaction from the metadynamics.py
6+ , the Walden inversion
7+ '''
8+
9+ import shutil
10+ import tempfile
11+ from pathlib import Path # a more Pythonic alternative to the os.path
12+ here = Path (__file__ ).parent
13+ # to the directory where the pseudopotential and orbital files are stored
14+ # In your case you change to the appropriate one
15+ pporb = here .parent .parent .parent / 'tests' / 'PP_ORB'
16+
17+ from ase .io import read
18+ from ase .atoms import Atoms
19+ from ase .optimize import BFGS , FIRE
20+ from ase .mep import NEB
21+ from ase .constraints import FixCartesian
22+ from abacuslite import AbacusProfile , Abacus
23+
24+ aprof = AbacusProfile (
25+ command = 'mpirun -np 16 abacus_2p' ,
26+ pseudo_dir = pporb ,
27+ orbital_dir = pporb ,
28+ omp_num_threads = 1 ,
29+ )
30+ PSEUDOPOTENTIALS = {
31+ 'C' : 'C_ONCV_PBE-1.0.upf' ,
32+ 'H' : 'H_ONCV_PBE-1.0.upf' ,
33+ 'F' : 'F_ONCV_PBE-1.0.upf' ,
34+ }
35+ BASISSETS = {
36+ 'C' : 'C_gga_8au_100Ry_2s2p1d.orb' ,
37+ 'H' : 'H_gga_8au_100Ry_2s1p.orb' ,
38+ 'F' : 'F_gga_7au_100Ry_2s2p1d.orb' ,
39+ }
40+ DFTPARAM = {
41+ 'calculation' : 'scf' ,
42+ 'nspin' : 2 , # because we are studying a reaction involving radical
43+ 'basis_type' : 'lcao' ,
44+ 'ks_solver' : 'genelpa' ,
45+ 'ecutwfc' : 60 , # good for example, wrong for production
46+ 'symmetry' : 1 ,
47+ 'gamma_only' : True ,
48+ 'init_chg' : 'auto'
49+ }
50+
51+ def impose_constraint (atoms : Atoms ) -> Atoms :
52+ '''imposing the constraint such that the C atom is fixed,
53+ H and F that leaves and arrives are fixed such that can only
54+ move along Z-axis'''
55+ # reset the constraint
56+ atoms .set_constraint (None )
57+ atoms .set_constraint ([FixCartesian (a = [0 ]),
58+ FixCartesian (a = [4 , 5 ], mask = (True , True , False ))])
59+ return atoms
60+
61+ def relax (atoms : Atoms ) -> Atoms :
62+ '''relax the structure to the minimum energy configuration'''
63+ # we do additional fully-constraint on the No.1 atom, and XY-constraint
64+ # on the last two atoms
65+ print ('Relaxation of atoms:' , atoms )
66+ with tempfile .TemporaryDirectory () as tmpdir :
67+ calc = Abacus (
68+ profile = aprof ,
69+ directory = tmpdir ,
70+ pseudopotentials = PSEUDOPOTENTIALS ,
71+ basissets = BASISSETS ,
72+ inp = DFTPARAM
73+ )
74+ atoms .calc = calc
75+ opt = BFGS (atoms , logfile = '-' )
76+ opt .run (fmax = 0.05 )
77+ print ('Done' )
78+ return atoms
79+
80+ # relax the initial and final states
81+ ch4f = ''' 6
82+ Lattice="15.0 0.0 0.0 0.0 15.0 0.0 0.0 0.0 15.0" Properties=species:S:1:pos:R:3
83+ C -0.00000000 -0.00000000 -1.50074532
84+ H 0.00000000 -1.01539888 -1.16330635
85+ H -0.87936122 0.50769944 -1.16330635
86+ H 0.87936122 0.50769944 -1.16330635
87+ H -0.00000000 -0.00000000 -2.57055225
88+ F 0.00000000 -0.00000000 0.88279136
89+ '''
90+ with tempfile .NamedTemporaryFile (mode = 'w' , suffix = '.extxyz' ) as f :
91+ f .write (ch4f )
92+ f .flush ()
93+ ini = relax (impose_constraint (read (f .name )))
94+
95+ ch3fh = ''' 6
96+ Lattice="15.0 0.0 0.0 0.0 15.0 0.0 0.0 0.0 15.0" Properties=species:S:1:pos:R:3
97+ C -0.00000000 -0.00000000 -1.50074532
98+ H 0.00000000 -1.01539888 -1.16330635
99+ H -0.87936122 0.50769944 -1.16330635
100+ H 0.87936122 0.50769944 -1.16330635
101+ H -0.00000000 -0.00000000 -3.57055225
102+ F 0.00000000 -0.00000000 -0.50000000
103+ '''
104+ with tempfile .NamedTemporaryFile (mode = 'w' , suffix = '.extxyz' ) as f :
105+ f .write (ch3fh )
106+ f .flush ()
107+ fin = relax (impose_constraint (read (f .name )))
108+
109+ # prepare images
110+ n_replica , images = 5 , []
111+
112+ for irep in range (n_replica ):
113+ calc = Abacus (
114+ profile = aprof ,
115+ directory = here / f'neb-idpp-replica-{ irep } ,
116+ pseudopotentials = PSEUDOPOTENTIALS ,
117+ basissets = BASISSETS ,
118+ inp = DFTPARAM
119+ )
120+ if irep == 0 :
121+ ini .calc = calc
122+ images .append (impose_constraint (ini ))
123+ elif irep == n_replica - 1 :
124+ fin .calc = calc
125+ images .append (impose_constraint (fin ))
126+ else :
127+ rep = ini .copy () if irep <= n_replica // 2 else fin .copy ()
128+ rep .calc = calc
129+ images .append (impose_constraint (rep ))
130+
131+ neb = NEB (images )
132+ neb .interpolate ('idpp' )
133+
134+ # optimize the band
135+ qn = FIRE (neb , trajectory = 'neb-idpp.traj' , logfile = '-' )
136+ qn .run (fmax = 0.05 )
0 commit comments