#!/bin/sh

# Example: Delta MEP-NOCV analysis of the donor-acceptor bond in ammonia borane, H3N-BH3
#
# The system is divided into two fragments, NH3 (electron donor) and BH3 (electron
# acceptor), both in the geometry they have in the complex. The complex is computed
# with the ETS-NOCV analysis and the DeltaMEP option, which fits the deformation
# density of every selected NOCV pair with ZlmFit and stores the fits on the ADF result file
# (ZlmFit_NOCV_<k> sections). The densf program then evaluates, for NOCV pair 3,
# the corresponding electrostatic potential (Delta MEP) on a grid,
# and writes it to a TAPE41 file that can
# be visualised in AMSview.
#
# Experimental geometry taken from
#   L. R. Thorne, R. D. Suenram, F. J. Lovas, J. Chem. Phys. 78, 167 (1983),
#   doi:10.1063/1.444528
#
#
# ---------------------------------------------------------------- fragment NH3
AMS_JOBNAME=nh3frag $AMSBIN/ams << eor
Task SinglePoint

System
   Atoms
      N   0.000000   0.000000   1.657600
      H   0.475571   0.823713   2.009061
      H  -0.951142   0.000000   2.009061
      H   0.475571  -0.823713   2.009061
   End
End

Engine ADF
   Basis
      Type TZP
      Core None
   End
   XC
      GGA BP86
   End
   Relativity
      Level None
   End
   Symmetry NOSYM
EndEngine
eor

# ---------------------------------------------------------------- fragment BH3
AMS_JOBNAME=bh3frag $AMSBIN/ams << eor
Task SinglePoint

System
   Atoms
      B   0.000000   0.000000   0.000000
      H   1.176251   0.000000  -0.308364
      H  -0.588126   1.018664  -0.308364
      H  -0.588126  -1.018664  -0.308364
   End
End

Engine ADF
   Basis
      Type TZP
      Core None
   End
   XC
      GGA BP86
   End
   Relativity
      Level None
   End
   Symmetry NOSYM
EndEngine
eor

# ---------------------------------------------------------------- complex NH3-BH3
# ETSNOCV computes the NOCV pairs, DELTAMEP fits the pair deformation densities.
# ENOCV 0.01 : only pairs with |eigenvalue| > 0.01 are analysed and fitted
#              (densf can use only the pairs that are stored on the result file).

AMS_JOBNAME=complex $AMSBIN/ams << eor
Task SinglePoint

System
   Atoms
      B   0.000000   0.000000   0.000000   adf.f=BH3
      N   0.000000   0.000000   1.657600   adf.f=NH3
      H   1.176251   0.000000  -0.308364   adf.f=BH3
      H  -0.588126   1.018664  -0.308364   adf.f=BH3
      H  -0.588126  -1.018664  -0.308364   adf.f=BH3
      H   0.475571   0.823713   2.009061   adf.f=NH3
      H  -0.951142   0.000000   2.009061   adf.f=NH3
      H   0.475571  -0.823713   2.009061   adf.f=NH3
   End
End

Engine ADF
   Basis
      Type TZP
      Core None
   End
   XC
      GGA BP86
   End
   Relativity
      Level None
   End
   Symmetry NOSYM

   Fragments
      BH3 bh3frag.results/adf.rkf
      NH3 nh3frag.results/adf.rkf
   End

   ETSNOCV
      Enocv 0.01
      DeltaMEP
   End
EndEngine
eor

# ---------------------------------------------------------------- densf
# Evaluate on a grid, for NOCV pair 3 only:
#   - the Coulomb potential of its fitted deformation density, Delta MEP.
# The pair must have been fitted by ADF; otherwise DENSF stops with an error.
# The field written to the output file is DeltaMEPNOCV%Coulpot_3 (restricted calculation).
# The deformation density itself can be requested separately with SUM 3 in the NOCV block.

AMS_JOBNAME=densf $AMSBIN/densf << eor
ADFFile complex.results/adf.rkf
OutputFile densf_dmep.t41
Grid Medium
End
NOCV
   DeltaMEP 3
End
eor

# Check that only the selected pair is written and that its grid values are finite.
"$AMSBIN/amspython" <<'PY'
from math import isfinite
from typing import List, Set

from scm.plams import KFFile

result: KFFile = KFFile("densf_dmep.t41")
fields: Set[str] = result.get_skeleton()["DeltaMEPNOCV"]
assert fields == {"Coulpot_3"}, fields
npoints: int = result.read_int("Grid", "total nr of points")
for field in sorted(fields):
    values: List[float] = result.read_reals("DeltaMEPNOCV", field)
    assert len(values) == npoints, field
    assert all(isfinite(value) for value in values), field
    assert any(value != 0.0 for value in values), field
PY