1from math import radians, sin, cos 2 3import pytest 4 5from ase import Atoms 6from ase.neb import NEB 7from ase.constraints import FixAtoms 8from ase.optimize import QuasiNewton, BFGS 9from ase.visualize import view 10 11 12@pytest.mark.calculator('nwchem') 13def test_h3o2m(factory): 14 # http://jcp.aip.org/resource/1/jcpsa6/v97/i10/p7507_s1 15 doo = 2.74 16 doht = 0.957 17 doh = 0.977 18 angle = radians(104.5) 19 initial = Atoms('HOHOH', 20 positions=[(-sin(angle) * doht, 0, cos(angle) * doht), 21 (0., 0., 0.), 22 (0., 0., doh), 23 (0., 0., doo), 24 (sin(angle) * doht, 0., doo - cos(angle) * doht)]) 25 if 0: 26 view(initial) 27 28 final = Atoms('HOHOH', 29 positions=[(- sin(angle) * doht, 0., cos(angle) * doht), 30 (0., 0., 0.), 31 (0., 0., doo - doh), 32 (0., 0., doo), 33 (sin(angle) * doht, 0., doo - cos(angle) * doht)]) 34 if 0: 35 view(final) 36 37 # Make band: 38 images = [initial.copy()] 39 for i in range(3): 40 images.append(initial.copy()) 41 images.append(final.copy()) 42 neb = NEB(images, climb=True) 43 44 def calculator(): 45 return factory.calc( 46 task='gradient', 47 theory='scf', 48 charge=-1 49 ) 50 51 # Set constraints and calculator: 52 constraint = FixAtoms(indices=[1, 3]) # fix OO 53 for image in images: 54 image.calc = calculator() 55 image.set_constraint(constraint) 56 57 # Relax initial and final states: 58 with QuasiNewton(images[0]) as dyn1: 59 dyn1.run(fmax=0.10) 60 with QuasiNewton(images[-1]) as dyn2: 61 dyn2.run(fmax=0.10) 62 63 # Interpolate positions between initial and final states: 64 neb.interpolate() 65 66 for image in images: 67 print(image.get_distance(1, 2), image.get_potential_energy()) 68 69 with BFGS(neb) as dyn: 70 # use better basis (e.g. aug-cc-pvdz) for NEB to converge 71 dyn.run(fmax=0.10) 72 73 for image in images: 74 print(image.get_distance(1, 2), image.get_potential_energy()) 75