Optimise and view the structure of ammonia dimer

[1]:

from ase.collections import s22

[6]:

name='Ammonia_dimer'
atoms = s22[name]

[7]:

type(atoms)

[7]:

ase.atoms.Atoms

[9]:

import ase.visualize
ase.visualize.view(atoms, viewer='x3d')

[9]:

ASE atomic visualization

[11]:

from pymolpro import ASEMolpro
from ase.optimize import BFGS
atoms.calc = ASEMolpro(ansatz='df-lmp2/aug-cc-pVTZ')
with BFGS(atoms) as opt:
    opt.run(fmax=0.0001)

      Step     Time          Energy          fmax
BFGS:    0 17:13:35    -3072.767830        0.039808
BFGS:    1 17:13:51    -3072.767945        0.031691
BFGS:    2 17:14:08    -3072.768093        0.038783
BFGS:    3 17:14:24    -3072.768165        0.024396
BFGS:    4 17:14:41    -3072.768200        0.014780
BFGS:    5 17:14:57    -3072.768213        0.012905
BFGS:    6 17:15:14    -3072.768237        0.015387
BFGS:    7 17:15:31    -3072.768293        0.029246
BFGS:    8 17:15:48    -3072.768409        0.044499
BFGS:    9 17:16:05    -3072.768583        0.048781
BFGS:   10 17:16:22    -3072.768728        0.031805
BFGS:   11 17:16:38    -3072.768776        0.011471
BFGS:   12 17:16:56    -3072.768786        0.007871
BFGS:   13 17:17:12    -3072.768791        0.009934
BFGS:   14 17:17:29    -3072.768798        0.010301
BFGS:   15 17:17:45    -3072.768809        0.008511
BFGS:   16 17:18:03    -3072.768821        0.009436
BFGS:   17 17:18:20    -3072.768835        0.009928
BFGS:   18 17:18:37    -3072.768852        0.012732
BFGS:   19 17:18:54    -3072.768871        0.012696
BFGS:   20 17:19:12    -3072.768885        0.007342
BFGS:   21 17:19:29    -3072.768890        0.001984
BFGS:   22 17:19:46    -3072.768890        0.000276
BFGS:   23 17:20:03    -3072.768890        0.000033

[13]:

def closest_interatomic(atoms):
    distances = atoms.get_all_distances()
    n = len(distances) // 2
    result = 1000.0
    for i in range(n):
        for j in range(n):
            result = min(distances[i][j + n], result)
    return result

[15]:

print("H-bond length:",closest_interatomic(atoms))

H-bond length: 2.1971031666750096

[30]:

from pymolpro import Project
p = Project(name)
import ase.io
ase.io.write(p.filename() + '/initial.xyz', s22[name])
p.write_input("""
geometry=initial.xyz
basis,aug-cc-pvtz
df-rhf; df-lmp2; optg,gradient=0.00001; put,xyz,final.xyz""")
p.run(wait=True)
assert p.status == 'completed'
method=p.xpath('//summary/@overall_method')[-1]

[32]:

final = ase.io.read(p.filename() + '/final.xyz')

[33]:

print("Initial H-bond length:",closest_interatomic(s22[name]),"\n"+method+" H-bond length:",closest_interatomic(final))

Initial H-bond length: 2.149250106516224
DF-LMP2/aug-cc-pVTZ//DF-LMP2/aug-cc-pVTZ H-bond length: 2.1971082552648555

[24]:

ase.visualize.view(final, viewer='x3d')

[24]:

ASE atomic visualization